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earth  was  developed  by  Biot  and  later  extended  by  Stoll  primarily  for 
describing  sound  propagation  in  ocean  sediments.  The  Biot-Stoll  model 
describes  sound  propagation  in  a  medium  composed  of  fluid  saturated  pores 
and  an  elastic  frame  (matrix).  The  differential  equations  which  arise  from 
this  formalism  allow  for  a  fast  wave  which  propagates,  primarily,  in  the 
frame  and  a  slow  wave  which  propagates,  primarily,  in  the  fluid.  For  our 
application  we  assume  that  the  fluid  is  air  and  that  the  porous  material  is 
a  layer  bewteen  the  air  (an  upper  semi-infinite  half  space)  and  homogenous 
clay  (a  lower  semi -inf inite  half  space). •>.  The  resultant  boundary  value 
problem  was  solved  for  amplitude  and  phase-  of  the  acoustic  wave  above  the 
surface,  acoustic  wave  in  the  pores  of  the  porous  layer,  seismic  motion  of 
the  pore  frame,  and  seismic  motion  of  the  lower  clay  half  space.  Shear 
motion  was  included  for  completeness.  The  incoming  sound  wave  was  assumed 
to  be  planar. 

'Calculations  suggest  that  the  predicted  transfer  functions  are  not 
strongly  affected  by  the  flow  resistance  of  the  porous  layer,  the  elastic 
properties  of  the  clay,  and  the  thickness  of  the  porous  layer.  These  three 
quantities  were  independently  measured./^ The  transfer  function  was  measured 
for  a  microphone  probe  buried  in  the  porous  layer,  and  triaxial  geophones 
in  the  porous  layer  and  underlying  clayJ  Measurements  were  made  for  acous¬ 
tic  frequencies  between  30  and  300  Hz*  for  angles  of  incidence  for  the 
acoustic  wave  between  5  and  20  ,  and  various  depths  for  buried  geophones 
and  microphone  probes.  j 

Experimental  values  of  the  various  transfer  functions  agree  with 
predictions  within  a  factor  of  two.  The  experimental  results  show  fine 
structure  not  predicted  which  we  attribute  to  multiple  layers.  Other 
quantities  predicted  by  the  model  (attenuation  of  signal  with  depth,  little 
dependence  on  angle  of  incidence,  acoustic  surface  impedance)  all  agree 
equally  well  with  experiment. 

We  conclude  that  the  Biot-Stoll  model  is  a  good  physical  description 
of  acoustic  to  seismic  energy  transfer.  Additional  avenues  for  future 
research  are  also  discussed.  Computer  programs  are  included  along  with  a 
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1.  INTRODUCTION 


The  objective  of  this  study  was  to  establish  the  physical  principles 

involved  in  coupling  airborne  sound  into  earth  motion.  The  study  was 

motivated  by  measurements  made  at  the  U.S.  Army  Corps  Engineers  Waterways 

Experiment  Station  (WES)  in  1977.  In  these  measurements,  a  speaker  was 

suspended  from  a  crane,  the  sound  level  in  the  far  field  was  measured  with 

a  microphone  above  the  surface,  and  geophone  response  was  measured  below 

1  2 

the  surface.  The  major  findings  *  were: 

1.  The  geophone  response  had  a  magnitude  1000  t;imes  larger  than  one 
would  expect  from  calculations  based  upon  simple  acoust  transmission 
through  a  boundary  between  two  perfect  fluids. 

2.  The  magnitude  of  the  geophone  response  is,  withii  ^  accuracy  of 
the  measurements,  independent  of  angle  of  incidence. 

3.  The  transit  time  from  speaker  to  microphone  was  approximately  the 
same  as  from  speaker  to  geophone  and  the  speed  of  propagation  was,  within 
the  accuracy  of  the  measurements,  the  same  as  the  speed  of  sound  in  air 
( 343  m/sec). 

4.  The  geophone  signal  varied  in  phase  as  one  would  expect  of  a 
Rayleigh  Wave. 

5.  The  geophone  signal  had  well  defined  maxima  and  minima  which  could 
be  accurately  explained  assuming  interference  in  the  upper  boundary  layer 
of  the  earth.  The  depth  of  this  layer  was  defined  by  independent  seismic 
refraction  studies. 

During  the  course  of  the  study  reported  on  here,  a  physical  model 
consistent  with  these  observations  has  been  developed  from  which  quantita¬ 
tive  predictions  can  be  made. 

The  first  inquiry  this  laboratory  received  concerning  this  problem  was 
from  WES  and  was  prompted  by  a  field  test  during  which  airborne  sound 


sources  had  detonated  smart  mines.  The  first  goal,  then,  was  to  develop  an 
algorithm  which  would  allow  the  sensor  to  discriminate  against  airborne 
sound.  This  was  done  by  discriminating  against  geophone  signals  which 
correlate  with  above  ground  microphone  signals.  More  recently,  this  corre¬ 
lation  has  been  used  to  identify  airborne  sound  source  in  a  security  system 
developed  by  Sandia  Laboratories.  The  coupling  between  airborne  and  seis¬ 
mic  signals,  then,  is  not  necessarily  detrimental  to  military  applications. 

In  the  following  sections  several  experiments  and  theoretical  results 
will  be  presented  which  represent  the  progression  of  events  leading  to  a 
physical  model  for  the  coupling  process. 

Initially,  the  theoretical  approach  was  to  treat  the  boundary  as  a 

separation  between  homogeneous  half  spaces  and  treat  the  shallow  porous 

345 

upper  layer  of  the  earth  as  an  impedance  matching  layer.  ’  ’  It  was 
thought  that  this  approach  plus  consideration  of  Rayleigh  Waves  in  the 
earth  instead  of  compressional  waves  would  give  the  measured  seismic  ampli¬ 
tudes.  Even  with  the  relatively  slow  Rayleigh  Wave  velocity,  however,  the 
acoustic  resistance  (pc)  ratio  between  air  and  the  earth  is  so  small  that 
most  all  the  impinging  energy  would  be  reflected.  This  approach  was 
abandoned  after  one  year. 

The  initial  experimental  concerns  were  the  extent  to  which  the  speaker 
was  decoupled  from  the  earth  and  the  potential  contributions  of  the  geo¬ 
phone  signal  by  direct  airborne  waves.  Several  experiments  were  conducted 
to  insure  direct  coupling  though  speaker  supports  would  not  account  for  the 
geophone  signal.  Speakers  were  isolated  from  the  ground  using  flexible 
cables  and  inner  tubes  with  no  effect  on  the  measured  seismic  signal. 
Perhaps  the  most  conclusive  experiment  was  the  firing  of  a  freely  falling 
blank  pistol  -  the  seismic/acoustic  coupling  remained  unchanged. 

At  this  stage  in  the  research  (about  one  year  after  beginning),  it  was 
concluded  that  airborne  sound  is,  indeed,  coupled  into  motion  of  the  earth 
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frame  with  an  amplitude  larger  than  can  be  predicted  from  homogeneous  layer 
theory. 

A  major  breakthrough  in  the  understanding  of  this  problem  came  as  a 
result  of  progress  on  another  research  project  in  this  laboratory.  To 
explain  measurements  of  sound  propagation  through  ocean  sediments,  it  is 
common  to  use  what  is  known  as  the  Biot-Stoll  Model. ^  This  model  basically 
allows  for  two  waves  in  the  porous  medium,  one  associated  with  the  matrix, 
another  associated  with  motion  of  the  pore  fluid  relative  to  the  matrix. 
At  a  meeting  organized  by  this  laboratory  on  the  status  of  long  range 
acoustic  propagation  research,^  the  possibility  of  using  the  Biot-Stoll 
Model  to  explain  coupling  of  airborne  sound  into  the  earth  was  explored. 
The  physical  principles  were  consistent  with  this  problem,  but  the  approxi¬ 
mations  used  in  underwater  propagation  codes  were  of  concern.  Following 
this  meeting.  Dr.  Attenborough  of  the  Open  University  of  England  agreed 
that  this  model  offered  most  promise  of  explaining  the  seismic/acoustic 
coupling  results.  It  was  agreed  that  a  joint  effort  would  be  fruitful,  and 
cooperative  measurements  in  this  laboratory  the  following  summer  (1982) 
were  arranged. 

During  the  time  Dr.  Attenborough  was  visiting  this  laboratory  it  was 
agreed  to  address  two  major  questions: 

1.  Is  there  an  acoustic  wave  in  the  earth  which  can  be  attributed  to 
motion  of  the  pore  fluid  relative  to  the  matrix? 


acoustic  wave  attributable  to  motion  of  air  in  the  pores  was  detected,  the 
attenuation  of  this  wave  with  depth  was  measured  and  favorable  comparisons 
with  calculations  made  assuming  a  rigid  porous  medium  were  shown.  These 
results  were  published  in  Reference  8.  The  same  model  used  to  compute 
attenuation  was  also  used  to  predict  the  velocity  of  the  acoustic  wave  in 
the  earth.  It  was  gratifying  to  find  that  the  predicted  speed  decreases 
rapidly  with  decreasing  frequency  such  that  for  a  typical  soil,  the  wave¬ 
length  becomes  independent  of  frequency  below  about  500  Hz  with  a 
wavelength  of  a  few  tens  of  centimeters.  Measurements  with  the  probe 
microphone  illustrated  that  in  the  same  medium,  an  acoustic  signal  could  be 
observed  down  to  depths  of  a  few  tens  of  centimeters.  The  pore  depth  is, 
in  fact,  comparable  to  the  acoustic  wavelength  even  at  very  low 
frequencies. 

These  results  reinforced  confidence  in  the  Biot-Stoll  Model  and  showed 
that  pervious  knowledge  in  acoustic  wave  propagation  in  ocean  sediments  and 
above  the  surface  of  the  earth  could  be  applied  directly  to  the 
seismic/acoustic  coupling  problem.  The  next  step  was  to  apply  the  Biot- 
Stoll  Model  quantitatively  dropping  the  earlier  assumption  (commonly  made 
in  acoustics)  that  the  frame  is  rigid.  The  result  of  that  effort  will  be 


presented  in  Chapter  2 


2.  THEORY 


In  this  section  application  of  the  Biot-Stoll  Model  of  sound  propaga¬ 
tion  in  a  porous  media  to  the  air/earth  boundary  will  be  described.  In  the 

first  section,  the  procedure  for  computing  propagation  coefficients  in  the 
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media  as  developed  by  Attenborough  *  will  be  presented.  In  the  next 
section  the  boundary  value  problem  for  this  case  will  set  up  and  solved. 
This  treatment  is  similar  to  that  presented  by  Sides ^  except  a  shear  wave 
term  is  also  included.  Finally  measurable  quantities  will  be  defined  and 
calculated. 

Absent  from  this  chapter  are  treatments  of  a  spherically  spreading 
source  (plane  waves  are  assumed),  impulsive  sources,  and  a  soil  with 
porosity  varying  with  depth.  The  effect  of  spherical  waves  could  be  to 
introduce  a  sub-surface  wave  similar  to  the  surface  wave  known  in  above 
ground  propagation  and  which,  above  the  ground,  is  dominant  for  certain 
source  receiver  geometries.  By  transforming  from  frequency  to  time  domain, 
the  current  theory  should  be  applicable  to  impulsive  sounds  but  that  pro¬ 
cess  has  not  yet  been  developed.  Varying  porosity  could  improve  coupling 
into  the  earth  frame  by  providing  a  gradual  change  in  the  acoustic 
resistance  with  depth.  These  two  extensions  to  the  theory  will  probably  be 
necessary  before  all  experimental  results  can  be  explained. 

The  Biot-Stoll  Model  for  wave  propagation  in  a  porous  elastic  medium 
was  used  to  describe  the  coupling  of  an  airborne  acoustic  wave  into  the 
porous  elastic  earth.  The  first  layer  or  few  tens  of  centimeters  of  the 
earth's  surface  is  porous  and  will  allow  for  the  propagation  of  an  acoustic 
wave,  whether  the  surface  is  sand  or  grass  covered  soil.  The  porosity  of  a 
sand  results  from  the  packing  of  the  sand  particles.  A  grass  covered 
soil's  porosity  might  result  from  decaying  organic  matter  and  grass  roots 
as  well  as  the  particle  packing.  This  first  layer  was  modeled  as  a  two 


dimensional  homogeneous  porous  elastic  medium  characterized  by  a  frame  or 
matrix  filled  with  air  overlying  a  semi-infinite  nonporous  elastic  clay. 


From  a  plane  wave  analysis  of  the  Biot-Stoll  differential  equations 
which  govern  the  propagation  of  compressional  waves,  the  attenuation  coef¬ 
ficients  and  phase  velocities  of  the  allowed  wave  motions  can  be 
determined.  After  applying  the  boundary  conditions  at  the  boundaries  of 
the  air-sand-clay  or  air-soil-clay  layer  to  solve  for  the  amplitudes  of  the 
matrix  and  fluid  motions,  the  specific  normal  acoustic  surface  impedance, 
the  attenuation  coefficients,  and  the  surface  transfer  functions  for  the 
porous  medium  will  be  found. 


Determination  of  Propagation  Constants 


The  Biot-Stoll  differential  equations  for  compressional  wave  propaga¬ 


tion  are 


V2(He-CO  -  (pe-pf £) 


(la) 


and 


V2(Ce-M£) 


(pfe-m£) 


.  Jlli. 

<  3t 


(lb) 


The  terms  H,  C  and  M  are  frequency  dependent  constants  that  characterize 
the  porous  elastic  matrix  and  will  be  discussed  later  in  this  section.  The 
quantity  e  represents  the  matrix  dilitation,  e  *  V*vJ  where  it  is  the  matrix 
displacement.  The  dilitation  of  the  fluid  relative  to  the  matrix  is  £  and 

£  *  V  •  w 

where 

$  -  n(it-u) 

is  the  relative  fluid  velocity;  U  is  the  fluid  displacement  and  fl  is  the 
porosity.  The  density  of  the  matrix  and  the  fluid  are  identified  with  the 


i 


v. 


V* 


V,' 


t? 


terms  p  and  p^  while  n  and  <  are  the  dynamic  fluid  viscosity  and  the 

permeability  of  the  matrix.  The  term  jr  allows  for  damping  through 

viscous  drag  as  the  fluid  and  matrix  move  relative  to  each  other.  The  term 

m  is  a  dimensionless  parameter  that  accounts  for  the  fact  that  not  all  the 

fluid  moves  in  the  direction  of  the  macroscopic  pressure  gradient  because 

the  pores  do  not  all  run  normal  to  the  surface.  This  term  is  related  to 

the  Biot  added  mass  term  and  allows  for  damping  due  to  inertial  drag 

between  the  fluid  and  matrix  as  they  move  relative  to  each  other.  This 

2  2 

quantity  is  expressed  as  q  p  where  q  ,  the  tortuosity,  is  a  constant 
greater  than  or  equal  to  one  and  has  the  particular  value  l/cos0  when  the 
pores  are  cylindrical  and  are  inclined  at  0  to  the  surface  normal. 

The  terms  H,  C,  and  M  are  complex  constants  that  characterize  the 
elastic  response  of  the  matrix.  These  constants  are  determined  from  the 
bulk  and  shear  moduli  of  the  solid  material,  the  bulk  modulus  of  the  fluid 
and  the  bulk  modulus  of  the  matrix.  The  relationships  are  as  follows 


H 


[D-Kb]2  % 


Kr~Kb 

D-Kb 


M 


Kri 

D-Kb 


(2a) 

(2b) 

(2c) 


K  [l+0(f£  -1)]. 

r  K.f 


(2d) 


In  the  above  expressions,  is  the  bulk  modulus  of  the  solid  material, 

is  the  bulk  modulus  of  the  fluid  and  y  and  K  are  the  shear  and  bulk  moduli 

b 

of  the  matrix.  The  elastic  constants  can  be  made  complex  through  a  complex 
compressibility  of  the  pore  fluid.  However,  at  this  stage  they  were 
assumed  to  be  real.  The  bulk  and  shear  moduli  were  determined  from  the 
expressions  for  compressional  and  shear  wave  speeds  in  an  elastic  solid, 
which  are 


»*.  »'  m* 

-■  -V  . 


7 


v  7?: vr; v  ej y<'.» rj 


'.,  •?  »'.J  >.p  ’.'‘•yv.w'/mj.m  atwj ,.v*.i?,.,»v."~ 
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*• 
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*£■- 
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% 

>  V; 
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•  r 
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L, 


V 


t  * 

Aj  * 


hft 

$ 

w 


§ 


lG/p) 


1/2 


(3) 


and 


[  (B+  -j  G)  /p  ] 1/2 


(4) 


In  these  expressions,  B  and  G  are  the  bulk  and  shear  moduli  and  p  is  the 
density  of  the  medium.  For  the  porous  soil,  it  was  assumed  B  was  the  frame 
bulk  modulus,  G  was  the  frame  shear  modulus  and  p  was  the  frame  density. 
In  the  non-porous  clay,  it  was  assumed  the  bulk  modulus  B  was  the  bulk 
modulus  of  the  solid  material,  and  G  was  the  shear  modulus  of  the  same. 
These  assumptions  imply  that  the  porous  medium  is  composed  of  the  same 
material  as  the  clay  but  that  it  has  been  made  porous  by  some  means,  say 
weathering  or  cultivation. 

In  equation  (lb),  Stoll6  replaces  the  viscosity  by  a  viscosity  correc¬ 
tion  factor,  nF(X) ,  where  X  is  a  dimensionless  quantity  that  is  related  to 
the  thickness  of  the  viscous  boundary  layer  at  the  pore  wall.  Biot1* 

developed  expressions  for  F(X)  for  cylindrical  and  parallel  sided  pores  in 
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terms  of  the  fluid  viscosity  and  pore  diameter.  Attenborough  suggested 
that  the  dimensionless  parameter  X  for  an  arbitrary  pore  shape  is 


x  m  [8pfq2tu] 1/2 
n 

/i" 


(5) 


where  —  is  a  shape  factor  ratio,  <p  is  the  flow  resistivity  and  u  is  the 
angular  frequency.  The  static  shape  factor,  n,  relates  X  for  an  arbitrary 
pore  cross  section  to  that  of  a  cylindrical  pore  cross  section  while  the 
pore  shape  factor,  s,  accounts  for  the  fact  that  the  pores  may  have  varying 
cross  sections  along  their  lengths.  The  ranges  of  n  and  s  as  suggested  by 
Attenborough  are  1 . 5>s>l>n>0. 5.  The  viscosity  correction  factor  F(X)  is  a 
complex  function  and  is  written  as 


i  fe 


r-V 


-  *vvvcvrvr. 


J_  >  /IXT[/i  X] 
4  1/1X1 


with  T[/iX]  =  [/i  X] /Jq  [/i  Xl ,  where  and  are  the  zero  and  first  order 


Bessel  functions  of  the  first  kind. 


The  Biot-Stoll  equations  of  motion  admit  two  plane  wave  solutions  and 


thus  two  propagation  constants  or  wave  speeds.  The  waves  are  referred  to 


as  the  type  one  and  type  two  waves  or  simply  as  fast  and  slow  waves.  For 


each  wave  there  is  both  fluid  and  solid  displacement  and  each  wave  is 


distinguished  by  its  propagation  coefficients.  In  sands  and  soils  the  fast 


wave  primarily  moves  in  the  matrix  or  solid  component  while  the  slow  wave 


moves  primarily  in  the  fluid.  The  exact  relationship  between  the  fluid  and 


solid  displacement  will  be  determined  in  a  later  paragraph.  The  slow  wave 


is  a  diffusion  wave  and  it  is  highly  attenuated  and  very  dispersive  as  it 


propagates.  The  fast  wave  is  analogous  to  the  compressional  wave  that 


propagates  in  a  solid  and  it  is  relatively  unattenuated  and  not  dispersive. 


From  a  plane  wave  analysis  of  equations  (la)  and  (lb)  the  frequency 


equation  determining  the  propagation  constants  can  be  found.  If  the  plane 


waves  e  *  Aexp[i(i.x-wt)  ]  and  5  *  Bexp[i(ix-ut) ]  are  chosen  as  solutions  to 


equations  la  and  lb,  the  following  expressions 


(Hi.2— po)2)A  +  (pfa)2-Ci2)B  -  0 


(Ci.2-pf uj2)A  +  (mwZ-M£Z-io)F(X)^)B  -  0 


2  w„2 


are  found  after  substitution.  The  determinant  of  the  coefficients  must  be 


zero  so  one  can  write 


|Hi2-pu)2 


2  2 1 
pfu>  -Cl 


I  2  2 

|Ci  -p^u) 


ma)2-M£2-iti)F(X)^| 


There  are  two  complex  roots  of  this  equation  from  which  both  the  attenua¬ 


tion  and  phase  velocities  for  the  fast  and  slow  waves  can  be  found.  The 


IMTV? 


r.VCV'.V. 


solution  to  this  expression  must  be  done  numerically  because  of  the  Bessel 
functions  involved,  and  a  routine  is  included  in  Appendix  A. 1. 

From  equation  (7a)  we  can  write  down  the  relation  between  the  fluid 
and  matrix  wave  amplitudes  for  each  wave  type  as 


i  A 


2  2 
(Hi,  -pa)  ) 

2  2 
(pfco  -C£  ) 


The  Biot-Stoll  differential  equations  which  govern  the  propagation  of 
shear  waves  are  not  included  here  but  are  found  in  Reference  6  along  with 
the  characteristic  frequency  equation  for  the  propagation  constants.  In 
this  work  it  was  assumed  that  a  shear  wave  is  non-dispersive ,  and  the  wave 
number  was  determined  from  the  measured  shear  wave  velocity. 

Now  that  the  propagation  coefficients  are  known  the  boundary  value 
problem  can  be  considered. 


Boundary  Value  Problem 


The  first  few  tens  of  centimeters  of  the  earth's  surface  was  modeled 
as  a  two-dimensional  porous  elastic  medium  of  depth  d^  overlying  a  semi¬ 
infinite  non-porous  elastic  medium.  The  intermediate  medium  was  either  a 
sand  or  soil  and  the  lower  medium  was  a  non-porous  clay.  The  air-soil  or 
air-sand  interface  was  assumed  to  be  a  free  surface  and  the  lower  surface 
was  assumed  to  be  in  welded  contact  with  an  impermeable  membrane  between 
the  two  media.  The  air  is  allowed  to  flow  across  the  upper  boundary,  while 
it  is  not  allowed  to  penetrate  the  elastic  clay  below  the  porous  soil. 

For  this  physical  system,  the  boundary  conditions  have  been  developed 
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by  Deresiewiez  and  Sfalak.  The  boundary  conditions  at  the  upper  inter¬ 
face  are  (1)  continuity  of  normal  particle  velocity,  (2)  continuity  of 
fluid  pressure,  (3)  continuity  of  total  normal  stress,  and  (4)  continuity 


of  tangential  stress.  At  the  lower  interface  the  boundary  conditions  are 
(5)  continuity  of  normal  frame  velocity,  (6)  continuity  of  tangential  frame 


10 


t 


C\  c,' 

K 


velocity,  (7)  continuity  of  fluid  velocity,  (8)  continuity  of  total  normal 
stress,  and  (9)  continuity  of  tangential  stress. 

The.  pressure  and  stress-strain  relations  for  a  porous  medium  in  the 
Biot-Stoll^  notation  are: 


T«x  ■ 

xr,  -  H.-2u(e2+ex)-C5 

t  »  He-2p(e  +e  )-C£ 
zz  x  y 

t  -  UY 
xy  z 

•T  -  UY 
yz  x 

T  -  uy 
zx  y 


‘f 

j- 

where  t .  , 

A 

ij 

terms  e  , 
x’ 

where  t  are  the  total  stresses,  and  is  the  pore  fluid  pressure.  The 


x  3x 

e  „ 

y  3y 

3uz 

A  a  — - 

z  3z 

^ 

+  ir> 
■  1/2<i7  +  %F> 


I! 


It  is  informative  to  evaluate  these  relations  in  a  medium  which  the 
porosity  is  either  one  or  zero.  If  the  porosity  of  the  medium  goes  to  one, 
0-1,  then  the  medium  is  a  fluid  and  the  relation  for  the  pressure  reduces 
to 


v  V 

j  V» 

J  '■> 


pf  -  M5 


because  e  -  V  *u  -  0  since  there  is  no  solid  component  present.  In  this 
last  expression  for  P^,  the  elastic  constant,  M,  reduces  to  the  bulk 


modulus  of  the  fluid.  In  the  lowest  medium,  where  the  porosity  is  zero, 

the  stress-strain  relations  reduce  to  that  of  an  elastic  solid. 

Specifically,  the  expression  for  x  reduces  to 

zz 


x  ■  He  -  2p(e  +e  ) 

ZZ  X  z 


since  there  is  no  fluid  component  present.  In  this  last  expression,  the 
elastic  frame  constant,  H,  reduces  to 


H*Kr+jU. 


In  the  following  M  will  be  set  equal  to  the  bulk  modulus  when  n  ■  1  and, 
for  a  non-porous  elastic  material,  Q  =  0,  H  will  be  evaluated  using 
equation  (25). 

Next  the  types  of  allowed  wave  motion  are  chosen  and  then  the  boundary 
conditions  are  applied  at  the  appropriate  interfaces.  In  the  air,  incident 
and  reflected  waves  are  allowed  having  wave  vectors  Z^  and  Z^,  with  1 2.^  |  » 
|i.r|.  Type  one  and  type  two  dilitational  waves  and  a  shear  wave  are 
transmitted  into  the  porous  medium  with  respective  wave  vectors,  Z^, 
and  Zy  At  the  lower  interface  the  three  waves  are  reflected  and  have  wave 
vectors  tj,  Z^  and  ~Z'^  with  | |  -  | A* |  ,  |&2|  -  |i.£|  and  -  |££|.  In 
the  clay  only  a  shear  wave  and  compressional  wave  are  transmitted  with  wave 
vectors  Z^  and  Z ,  respectively.  Each  allowed  wave  is  associated  with  an 
angle  of  incidence,  reflection,  etc.  as  indicated  in  Figure  2.1.  In  the 
air  the  incident  and  reflected  waves  are 


*  ^B^expUCZ^cosS^x+J^sinO^z-uit)] 


U  ■  Z  B  exp[i(£  cos0  x-Z  sin0  z-wt)] 


• 


where  and  Bf  are  the  fluid  displacement  amplitudes.  In  the  porous 


medium, 

the  matrix  corapressional  waves  are 

*1 

=*  i^A^expf  i(  J^sin0  ^x+Sl^cos0  ^z-iut )  ] 

(28) 

*1 

CM 

iO 

■  £2A2exp[i(Jl2sine2x+Jl2cos02z-ut)  ] 

(29) 

f, 

•V 

•  a 

si 

=  £|A|exp[i()l|sin02x-)i^cos0  J^z-ojt)  ] 

(30) 

• 

•  j 

S2 

-  £2A2exp[i(£2sin02x-2.2cos02z-ut)  ] 

(31) 

where  u^  and  ^  are  the  type  one  and  two  waves  transmitted  into  the  porous 
medium  and  uj  and  u^  are  the  type  one  and  two  waves  reflected  at  the  lower 
boundary.  The  are  the  matrix  wave  displacement  amplitudes.  Associated 
with  each  matrix  wave  in  the  porous  medium  is  a  relative  fluid  wave  which 


is  expressed  as 

“  £  jB^expUCJl^sinQ  ^x+£^cos0 ^z-ut)]  (32) 

S2  *  £2B2exp^1^2sin92x+il2COS02Z_<i)t^  (33) 

w|  »  £|B|exp[i(&^sin0jLx-S.|cos9jZ-*wt)]  (34) 

■  i'B'expEKA'sine'x-il'cose'z-not)]  (35) 


where  the  B^  are  the  relative  fluid  wave  amplitudes.  In  the  porous  medium 
there  are  also  shear  waves  transmitted  and  reflected  which  are  expressed 
as 


*  8Eexp[i( ^sinO^x+Ji^cosO^z-wt) ]  (36) 

u^  *  S'E'exptKJl^sinO^x-J^cosB^z-wt)  ]  (37) 

where  E  and  E'  are  matrix  displacement  amplitudes.  The  transmitted  com- 
pressional  and  shear  waves  in  the  clay  are 

«  £^Atexp[i(JljSin0gX+£gCos0^z-a>t)  ]  (38) 

u®  ■  8"Etexp(i(2.^sin0^x+il^cos0^z-u)t)  ] 


(39) 


and  A(  and  are  again  the  displacement  amplitudes.  Using  these  plane 
waves,  the  boundary  conditions  at  the  two  interfaces  can  be  evaluated. 


Br  - 


Figure  2.1.  The  layer  model  used  indicating  notation  for  angles  of  inci¬ 
dence,  reflection  and  refraction,  wave  amplitudes,  propagation  vectors  and 
polarizations. 

At  the  upper  boundary,  the  first  boundary  condition  is  continuity  of 
normal  particle  velocity  and  is  expressed  as 

-  (l-«)u^+  (40) 

where  a  bar  is  used  to  indicate  wave  motion  in  the  upper  medium.  Using 
w  -  n(tr-U),  Uj^  can  be  eliminated  from  the  boundary  condition  to  yield 

Uj^  ■  Ujj-Wj^.  (41) 

Substituting  in  the  appropriate  waves  yields  the  following  equation 


U^cos0^-Urcos0 f  *  u^cos0 ^+U2co8®2“uicos®i 
-  u^cosO^-w^cos©  ^-w^cos©  2 
+  wjcos0  J+w^cos©  j+UjSin©  ^-u^sin©^. 


After  substituting  for  z  =  0,  Snell's  law  is  verified  yielding: 


This  boundary  condition  equation  can  be  rewritten  as 


(Bi~Br)cos0i  *  Aj ( 1-m^ )cos0  ^  +  A£( l-n^cosO ^ 

-  A|(l-m^)cos0 1  -  A2(l-m2)cos02 

+  Esin©^  -E'sinS^  (47) 

where  the  relative  wave  amplitudes  have  been  replaced  by  the  matrix  wave 
amplitudes,  using  m^  »  A^/B^.  The  second  boundary  condition  is  continuity 
of  fluid  pressure.  The  expression  for  the  pressure  is 


Pf  -  M£  -  Ce 


which  when  evaluated  at  the  boundary  yields 


VitBi+Br]  *  A1^1<C“®1M)  +  A2Z2(C-m2M)  +  AjJ^ (C-n^M) 

+  A^£2(c-m2M)* 


The  continuity  of  total  normal  stress,  xzz»  is  the  third  boundary  condition 
and  upon  substitution  yields  the  following  at  z  ■  0, 

Kf  ^i (Bi+Br)  ■  A^ 2,^  [H— m^C-2psin^0  ^  ]  +  A2£2  [H-n^C^^sin^l 
+  A|  [H-m^C-2ysin^ 0  ]  +  A^j ^ 

+  2EU&2sin®3cos93  +  2E 'M£2sin®3cos03*  (50) 

In  the  fluid,  t  was  set  equal  to  the  negative  of  the  pressure.  The  last 
zz 

boundary  condition  to  be  evaluated  at  the  upper  interface  is  continuity  of 
tangential  stress, 


(51) 


i  3ux  3u7 

t2*  ■  b  (ir +  ir> 


The  fluid  is  assumed  not  to  support  any  shears  and  thus  in  the  fluid 

was  set  equal  to  zero.  This  boundary  condition  results  in 


0  *  2A^£^sin0  ^cos0  ^  +  2A2^2s^-n® 2C0S® 2 

-  2Aj£^sin0  ^cos0  ^  -  2Aj£2sin02cos0 2 

-  EJl^cos^^-sin^^)  +  E'S,2(cos^03_sin^02). 


Next,  the  five  boundary  conditions  at  the  lower  interface  are  evaluated. 
The  fifth  boundary  condition  is  continuity  of  normal  frame  velocity, 


U1  ‘  ul* 


Substitution  of  the  plane  waves  at  the  lower  interface  yields 


u^cos0  j  +  u2cos02  -  u jcos0  ^  -  u2cos02  +  u^sin0^ 

c  s 

-  u'sinO.  »  u  cos0c + u  sin0, 

3  3  t  5  t  4 


After  setting  z  *  d  in  this  expression,  Snell's  law  is  again  verified 
equation  (54)  may  be  written  as 


A^cosO^exptiJ^cosO^]  +  A2cos02exp[i£2cos02d  ] 

-  AJcos0  jexpf-iJljCOsO^d]  -  A^cosO^xpUJ^cosO^d] 

+  EsinO^expfii^cosS^d]  -  E'sin03exp[-i£3cos03d] 
*  A^cos6^exp[i£^cos0^d]  +  Etsin0^exp[i£^sin0^d]. 


The  sixth  boundary  condition  is  continuity  of  tangential  frame  velocity, 


ut  -  dt. 


This  boundary  condition  yields 


u^sinO 1  +  u2sin02  +  u|sin02  +  u2sin@2  -  u^cos03  -  u3cose3 

c  s 

-  u  sin9c  -  u  cos0. 
t  5  t  4 
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and  at  z  =  d  reduces  to 


AjSin0  ^exp[i£^cos0  ^d]  +  A2sin02exp[i2,2cos0  2d] 

+  A|sin0 ^exp[-i£^cos0 ^d]  +  AjSinO 2exp[-i2,2C0S9 2^ 1 

-  Ecose^exptiJl^cos^d]  ”  E ’cosO^expf-iJl^cosO^d] 

*  Atsin0 ^exp(ii,gCos0gd]  -  Etcos0^exp[i£^cos8^d] .  (58) 

The  seventh  boundary  condition  is  continuity  of  normal  fluid  velocity, 
however,  since  the  fluid  cannot  penetrate  the  lower  medium,  the  normal 
fluid  velocity  must  go  to  zero  at  the  boundary,  w^  *  0.  Substitution  in 
this  expression  yields 

w^cos0^+w2cos02-w|cos0^-w^cos02  *  0.  (59) 

Using  the  relation  between  the  relative  wave  amplitude  and  the  matrix  wave 
amplitude,  equation  (59)  can  be  written  as,  at  z  •  d, 

A^m^cosO^exptiJl^cosO^d]  +  A2m2Cos0  2exp[i£.2cos92c^ 

-  Ajm^cosO  ^exp[-iJijCOS0  jd]  -  A2m2COs62exp[-i£2cos02<*J  *  0.  (60) 

Continuity  of  normal  stress,  Tzz»  is  the  eighth  boundary  condition  and  it 
yields,  at  z  -  d, 

2 

A^£^ [H-m^C-^psin  0  ^ ]exp(i£^cos0  ^d) 

2 

+  A2£2tH-m2C-2lIsin  92^exp^i)!'2COS02d^ 

+  Aj^  [H-m^C-2ysin29  ]exp(-i£  jCosO  ^d) 

+  A^^  [H-m^C-2y  sin20  ^  ]  exp(-i£  2cos9  2^ 

+  2p£2Esin93cos93exP(i^3cos93<i) 

+  2y£  ^E  *  sin©  ^cos©  ^exp^iX,^039  3d) 

2 

*  H’At£jexp[i£ jCos0 jd]  -  2p'At£jCos  0 ^exp[ i£^cos0 ^d ] 

-  2y  'E^^sinO^cosO^expUZ^cosO^d] ,  (61) 
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where  H*  and  p '  are  elastic  constants  in  the  clay.  The  last  boundary 


condition  is  continuity  of  tangential  stress,  t zx 
this  expression  yields,  at  z**d, 


x 


zx 


Substitution  in 


2pAj  Jl^sinQ^cose^expCi  ^cosO^d) 

+  2pA2^2si-n®2cos®2ex^^1';i2COS^2^^ 

-  2pAj  Jl^sinS^cosS^expC-iJl^cose^d) 

-  2uA^  i.2s^-n® 2C0S  ®2eX^^— ^  5'2COS  ® 2^  ^ 

2  2 

-  Ep£^(cos  0^-sin  0^ )exp(i 2.2C0S®3<i) 

2  2 

+  E'u£^(cos  0^-sin  ©^expC-ii^cosO^d) 

■  2Atp'  X.^sin0^cos0^exp(i£,gCos0^d) 

2  2 

-  E  u'Jl^Ccos  0^-sin  8^)exp(i£^cos0^d).  (62) 

These  nine  boundary  condition  equations  can  be  solved  simultaneously  to 
determine  the  amplitudes  of  the  allowed  types  of  wave  motion  relative  to 
the  incident  wave  amplitude.  From  these  amplitudes,  the  acoustic  propaga¬ 
tion  coefficients  of  the  porous  medium  as  well  as  the  surface  impedance  and 
transfer  functions  for  the  upper  boundary  can  be  found.  The  routine, 
Ampl. Fortran,  used  to  determine  those  amplitudes  calls  an  IMSL  subroutine. 
The  main  routine  is  included  in  Appendix  A. 2. 


Measurable  Quantities 


Predictions  of  seismic/acoustic  coupling  in  the  field  should,  when 

possible,  rely  only  upon  quantities  measurable  in  the  field.  Measurements 
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of  jacketed  and  unjacketed  bulk  moduli  are  interesting  in  the  laboratory 
but  not  practical  in  the  field.  One  goal  of  this  study  was  to  determine 
which  quantities  must  be  measured  in  order  to  predict  seimic/acoustic 
coupling. 

Acoustic  measurements  outdoors  usually  require  that  the  acoustic  impe¬ 
dance  of  the  surface  be  known.  Several  techniques  have  been  developed  for 


this  purpose 


It  has  been  found,  that  for  many  surfaces,  the  acoustic 


impedance  can  be  computed  from  the  flow  resistance.  The  flow  resistance 
can  be  readily  measured  in  the  field  or  estimated  by  an  on-site  evaluation 
of  the  surface.  During  the  course  of  these  measurements,  the  ability  to 
measure  the  attenuation  of  the  acoustic  wave  with  depth  in  the  soil  using  a 
buried  microphone  was  developed.  The  acoustic  attenuation  in  the  soil  is, 
typically,  easier  to  measure  than  the  acoustic  impedance  yet,  as  will  be 
shown  later,  carries  the  same  information  about  the  surface.  The  response 
of  a  buried  geophone  to  an  airborne  acoustic  wave  and  to  seismic 
disturbances  can  be  measured  in  the  field.  By  comparing  the  geophone  or 
probe  microphone  signal  to  the  signal  from  a  microphone  above  the  surface, 


one  can  determine  the  transfer  function  for  the  boundary.  Finally,  by 
mechanically  disturbing  the  surface,  the  seismic  velocities  can  be  measured 
as  well  as  layer  depths. 

Of  those  quantities  underlined  above,  flow  resistance,  seismic  veloci¬ 


ties,  and  layer  depths  must  be  experimentally  measured.  Experimental 
measurements  of  these  quantities  will  be  described  in  the  next  chapter. 
The  other  quantities  can,  in  principle,  be  computed  from  flow  resistance, 
seismic  velocities,  layer  depths,  and  tabulated  physical  properties  of  the 
soil.  In  the  following,  a  description  of  the  calculations  of  the  other 
underlined  quantities  will  be  given  using  the  Biot-Stoll  Model. 


Acoustic  Impedance  of  the  Surface 


The  normal  specific  acoustic  impedance  of  a  surface  is  the  ratio  of 
the  acoustic  pressure  to  the  particle  velocity  at  the  surface  divided  by 
the  acoustic  resistance  of  air,  Pqc0*  This  is  a  complex  quantity  and  is 
expressed  as 


z  ■  z  .  +  z. 

real  imag 


The  surface  pressure  in  terms  of  the  incident  and  reflected  wave  amplitudes 


P  ,  (to)  =  -iK,£J  [B.+B  ] 
surface  f  il  i  r 


and  the  velocity  at  the  surface  is 


VUurfa<!<!><“)  '  -*“cose1tBrBr'- 


The  normal  surface  impedance  is  then 


_  JP  Kf  l±  [Br+Bj] 
Z  "  V  “  aj()cos0i  [Br-Bi] 


_z _ Kf  [Bt+Br] 

zair  PqCq2cos0^  [Bj[-Br]  * 

In  Appendix  A. 3  is  a  subroutine,  Sufimp. Fortran,  used  to  calculate  the 
surface  impedance. 


Acoustic  Attenuation  in  the  Soil 


The  acoustic  attenuation  coefficients  of  the  porous  medium  are  found 
by  determining  the  magnitude  of  the  acoustic  pressure  at  two  different 
depths  below  the  surface.  The  acoustic  pressure  below  the  surface  is  the 
sum  of  the  fluid  pressure  due  to  the  fast  and  slow  waves.  Using  the  fast 
wave  amplitude  and  equation  (9),  for  example,  the  amplitude  of  the  fast 
relative  fluid  displacement  at  a  depth  z  is 


P-  .  -  iA. I , (C-m. M)exp( i£ . cos0 , z). 
last  11  1  1  1 


Expression  for  the  fluid  pressure  due  to  the  other  waves  in  the  porous 
medium  are 
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Pglow  *  iA2^2(C-m2M)exp(i£2cos92z) 

Pi  ■  iA! £ . (C-m, M)exp(-i£ , cos0  .  z) 
fast  11  1  r  1  1 

Pslow  =  iAj*.2(C-m2M)exp(-i£2cos02z) 


r;  »v  r./.'  r-  *\  vr,  r.  ir.  r.  v.  v.  ^  ^ 

(69) 

(70) 

(71) 


where  again  the  prime  indicates  a  wave  reflected  from  the  lower  boundary. 
The  total  fluid  pressure  at  a  depth  z_  below  the  surface  is  the  sum  of  the 
pressure  due  to  the  relative  fluid  waves  and  is  expressed  as 


(72) 


(73) 


P.  ,  (u>,z)  =  P,  +  P  ,  +  P'  +  P\  . 

below  fast  slow  fast  slow 


The  attenuation  coefficient  in  dB  per  unit  depth  is  then 

a  **  [£n|P  1  (w,z)|  -  £n|P,  1  (oo,z-)|]* 

1  below  *  '  1  below  2  1  (z^-z^) 

This  calculation  is  done  in  the  routine  Ampl. Fortran. 


Transfer  Function  of  the  Boundary 

In  this  section  the  acoustic  and  seismic  transfer  function  for  the 
air-soil  interface  will  be  determined.  The  acoustic  transfer  function  is 
the  ratio  of  the  acoustic  pressures  below  and  above  the  surface.  This  is  a 
complex  quantity  and  expressions  for  both  the  magnitude  and  phase  will  be 
included.  The  acoustic  pressure  at  a  height  z^  above  the  surface  is  due  to 
the  incident  and  reflected  wave  amplitudes  and  is  expressed  as 


P  ,  (w,z)  *  P.  (u),z)  +  P  .  ,(u>,z) 

above  incident  reflected 


(74) 


where 


Plncident(“’2)  ‘  -lKf*tBtex',(Ulc099lz) 


(75) 


and 


P  ,(ui,z)  ■  -iK,£  B  exp(-i£  cos0  z). 

reflected  f  r  r  *  r  r 


(76) 
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Using  the  expression  for  the  acoustic  pressure  below  the  surface,  equation 
(72),  the  acoustic  transfer  function  is  expressed  as 


,,  pbeW“>z) 
T(<i>  ,z ,z  )  =  Pabove(t0>z.) 


The  seismic  transfer  function  is  the  ratio  of  the  seismic  velocity 
response  of  the  matrix  below  the  boundary  to  the  acoustic  pressure  above 
the  boundary.  The  total  velocity  response  of  the  matrix  will  be  the  vector 
sum  of  the  frame  velocity  of  each  of  the  allowed  matrix  waves.  However, 
the  total  normal  velocity  response  of  the  matrix  will  be  determined  so  that 
it  can  be  later  compared  to  experimental  results.  The  total  normal 
velocity,  V  ,  at  a  depth  z  below  the  surface  is 

V„  •  V  +  a21  *  4il  +  S21  +  431  +  % 

V  -  -iuA,  cos8.  expCiJl,  cos6 ,  z)  +  itoA.'cosQ .  exp(-i£  cos0  z) 
n  11  11  ll  li 

-  iwA2cos02exp(iJl2cos02z)  +  io)A2cos02exp(-iil2cos02z) 

+  iwEsinO^expCiil^cosO^z)  +  iujE'sinOjexpC-iJ^cosO^z).  (78) 

An  expression  for  the  seismic  transfer  function  is 


T(o)  ,z  ,z '  ) 


Vn(oo.z) 

P above (u » z ' ) 


The  radial  transfer  function  as  well  as  an  air-clay  transfer  function  is  a 
simple  extension  of  this  work.  Routines  for  calculating  these  transfer 
functions  are  included  in  Appendix  A. 4,  A. 5,  and  A. 6. 
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3.  EXPERIMENTAL  MEASUREMENTS 
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During  the  course  of  this  study,  a  variety  of  experimental  measure¬ 
ments  were  performed  to  define  the  mechanism  responsible  for 
seismic/acoustic  coupling  and  insure  that  this  quantity  is  measured.  The 
individual  experiments  along  with  data  analysis  procedures  will  be 
described  in  the  following.  In  many  cases,  the  experiments  were  designed 
to  answer  specific  experimental  questions.  In  those  cases  the  results  and 
conclusions  will  be  presented  with  the  experimental  procedure.  The  major 
measurements  were  of  the  transfer  function.  Interpretation  of  those 
results  will  be  described  Chapter  5. 

Preliminary  Measurements 
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Considerable  effort  was  expended  to  make  certain  the  sound  system  was 
mechanically  decoupled  from  the  surface  of  the  earth  to  eliminate  direct 
mechanical  transfer  of  energy.  The  speaker  systems  were  either  suspended 
from  a  large  A  frame  using  various  elastic  straps,  or  supported  on  the  bed 
of  a  3/4-ton  truck  which  was  placed  on  inflated  rubber  bags. 

In  order  to  measure  the  degree  of  isolation,  two  types  of  tests  were 
conducted.  Impulses  on  the  side  of  the  speaker  cabinet  were  generated  and 
the  resulting  acoustic  and  seismic  waves  measured  using  a  collocated  geo¬ 
phone  and  microphone.  As  isolation  systems  were  applied  the  magnitude  of 
the  seismic  signal  was  observed  on  the  geophone.  The  microphone  provided  a 
time  scale  to  permit  separation  of  a  seismic  wave  generated  mechanically  at 
the  source  and  an  acoustically  coupled  wave.  In  all  cases  it  was  possible 
to  reduce  the  mechanical  coupling  below  the  ambient  background  seismic 
noise  level.  Geophone  response  to  sweep  tones  generated  by  the  speaker 
were  also  measured  as  the  decoupling  was  applied.  The  effectiveness  of  the 
decoupling  as  a  function  of  frequency  was  observed  in  this  way.  It  was 
possible  to  ascertain  that  no  specific  resonances  existed  in  the  isolation 
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system  which  would  selectively  pass  low  frequency  energy  to  the  earth  and 
geophone. 

A  geophone  isolated  from  the  ground  and  suspended  in  air  would  respond 
to  airborne  acoustic  energy.  Measurements  were  made  to  determine  if  acous¬ 
tic  waves  could  generate  a  geophone  response  when  the  geophone  was  placed 
in  a  variety  of  porous  materials,  but  still  isolated  from  the  earth  and 
seismic  energy  sources.  The  airborne  sound  produced  a  geophone  response 
when  the  geophone  was  covered  in  isolated  containers  of  plastic  beads,  of 
sands  and  gravels  with  a  variety  of  porosities.  The  overall  response  of 
the  geophone  varied  with  the  surrounding  materials,  but  in  each  case  the 
results  could  be  explained  in  terms  of  a  vibrating  mass  loaded  system 
driven  by  the  airborne  acoustic  energy.  The  data  indicated  no  measureable 
contribution  to  geophone  response  due  to  motion  of  the  geophone  relative  to 
the  matrix  material  when  driven  in  this  isolated  system. 

Measurements  of  the  Acoustic  Wave  in  the  Earth 

The  measurements  of  the  sound  pressure  below  the  surface  of  the  earth 
were  made  with  a  specially  designed  probe  as  shown  in  Figure  3.1.  It 
consists  of  an  AKG  microphone  and  preamplifier  housed  in  a  cylindrical 
brass  tube.  The  inner  cylinders  shown  were  used  to  increase  the  mass  of 
the  probe  and  to  seal  any  possible  air  leaks.  The  nose  cone  contains  small 
holes  to  allow  the  active  element  of  the  microphone  to  sense  the  acoustic 
field.  A  rubber  membrane  was  stretched  over  the  microphone  elements  to 
protect  it  from  any  particles  that  penetrate  the  nose  cone  when  the  probe 
was  inserted  in  a  predrilled  hole.  A  comparison  of  the  probe  microphone 
response  to  a  standard  AKG  microphone  is  shown  in  Figure  3.2.  The 
frequency  response  for  the  probe  is  essentially  flat  to  3kHz  where  it 
begins  to  roll  off  to  a  new,  less  sensitive,  plateau  near  7kHz.  The 
frequency  response  of  the  probe  is  essentially  unchanged  when  the  probe  is 


Brass  Cylinder 


Air  Holes 
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Figure  3.1.  Specially  designed  probe  microphone  used  to  measure  the  acous¬ 
tic  sound  field  in  the  porous  soil. 


FREQUENCY  (KHz) 


Figure  3.2.  A  comparison  of  probe  microphone  response  ( - )  to  a  standard 

AKG  microphone  ( - ).  Also  shown  is  the  response  of  the  probe  when  pointed 

away  from  the  sound  source. 


pointed  away  from  the  sound  source,  the  normal  geometry  for  measurements 
below  the  earth's  surface. 

Before  the  probe  was  inserted  into  the  ground,  a  drill  was  used  to 
bore  a  hole  1  inch  in  diameter  to  the  desired  depth.  The  probe  was  then 
forced  into  the  hole  and  STP  Oil  Treatment  was  poured  around  the  probe  near 
ground  level  to  further  seal  the  probe  entry.  This  procedure  worked  well 
giving  reproducible  results  for  sandy  soil  and  loam.  For  clays,  the  soil 
in  the  vicinity  of  the  probe  would  quite  often  crack  providing  air  paths 
not  previously  available.  This  cracking  gave  results  which  differed  from 
one  attempt  to  the  next.  It  was  generally  assumed  that  the  measurement 
giving  the  smallest  sound  pressure  at  the  probe  was  best. 

Occasionally,  the  holes  in  the  tip  of  the  probe  became  filled  with 
dirt  when  the  probe  was  inserted.  This  occurance  could  only  be  detected 
when  the  probe  was  removed  and  examined,  so  several  data  runs  had  to  be 
rejected  after  the  fact. 

The  results  of  an  extensive  series  of  probe  microphone  measurements 
are  published  in  Reference  8.  In  subsequent  measurements,  the  probe  micro¬ 
phone  output  was  recorded  along  with  geophone  response,  but  the  microphone 
depth  was  not  varied.  Referring  again  to  Figure  3.1,  when  the  probe  is 
inserted  in  the  ground  where  the  speed  of  sound  is  less  than  the  speed  in 
the  air,  the  probe  microphone  response  with  frequency  could  be  expected  to 
differ.  There  is  no  clear  theoretical  relation  to  predict  the  microphone 
response  when  it  is  inserted  into  the  earth  since  the  surrounding  medium  is 
changed  dramatically.  One  would  expect,  as  a  minimum,  that  the  frequency 
at  which  the  response  begins  to  roll  off  should  be  related  to  the  wave¬ 
length  of  sound  in  the  pores  which  scales  with  wavelength.  If  the  speed  of 
sound  in  the  pores  is  10  m/sec,  one  would  expect  the  microphone  response 
to  begin  to  roll  off  at  about  100  Hz.  This  feature  of  the  data  will  be 


demonstrated  below 


The  speed  of  sound  in  the  pores  was  measured  in  a  hemispherical  dish 
formed  in  the  ground  with  a  radius  of  0.5  m  and  filled  with  sand  with  a 
measured  flow  resistivity  of  85  rayls/cm.  Figure  3.3  shows  the  relative 


Figure  3.3.  Relative  positions  of  the  sound  source  and  sensors  used  to 
measure  acoustic  phase  velocity  and  attenuation. 

positions  of  the  source  and  sensors  used  in  the  experiment.  A  reference 
AK.G  microphone  was  positioned  approximately  5  cm  above  the  sand  and  the 
probe  microphone  inserted  to  an  approximate  depth  of  five  cm.  A  50  Hz 
acoustic  tone  was  broadcast  at  an  angle  of  incidence  of  approximately  20° 
from  the  horizontal.  The  relative  delay  between  the  probe  and  reference 
microphone  was  recorded  using  a  Nicolet  3091  digital  oscilloscope.  Next, 
the  probe  was  pushed  2.0  cm  further  down  in  the  sand  and  the  relative  delay 
again  recorded.  The  sound  speed  in  the  pores  is  the  ratio  of  the  change  in 
the  probe  depth  to  the  change  in  the  phase.  The  probe  was  pushed  down  in 
2.0  cm  increments  to  a  depth  of  11  cm  and  the  phase  difference  recorded  at 
each  probe  position.  This  experiment  was  repeated  at  several  frequencies 
between  50  Hz  and  500  Hz.  Results  of  the  average  phase  speed  in  the  pores 


are  plotted  in  Fig.  3.4.  The  attenuation  of  the  acoustic  wave  was  measured 
in  the  same  sand  using  sweep  tones,  results  are  shown  in  Figure  3.5.  The 
solid  line  was  calculated  using  equation  (73). 


Measurements  of  the  Transfer  Function 


Experimental  Configuration 


The  experimental  data  were  collected  at  two  different  locations,  the 
UM  test  field  and  Sardis  beach.  The  UM  test  field  is  two  miles  south  of 
Oxford,  MS  and  a  map  showing  the  field  size,  relative  elevations  and  tree 
lines  is  included  in  appendix  C.l.  Sardis  beach  is  located  at  John  W.  Kyle 
state  park  which  is  located  approximately  ten  miles  west  of  Oxford  and  a 
map  of  the  park  is  included  in  appendix  C.2.  Several  experiments  were 
performed  at  each  site  in  order  to  characterize  the  physical  properties 
necessary  for  theoretical  calculations. 

The  flow  resistivity  was  measured  using  an  apparatus  described  by 
Leonard^.  A  cylindrical  sample  of  earth  was  taken  from  the  ground  and  the 
rate  of  air  flow  through  the  sample  was  measured  as  a  function  of  differen¬ 
tial  pressure  across  the  sample.  Measurements  of  flow  resistance  were  made 
of  each  soil  sample  using  at  least  three  differential  pressures  across  the 
sample  .  Consistency  of  these  results  were  always  within  5%.  However,  it 
was  found  that  there  were  inconsistencies  in  the  measured  flow  resistances 
from  multiple  samples  from  the  same  site.  Values  of  flow  resistance  varied 
as  much  as  a  factor  of  two  at  the  UM  test  field,  however,  in  sand  these 
measurements  were  always  within  10%.  It  was  assumed  that  these  discrepan¬ 
cies  occured  because  of  variations  in  soil  structure  in  the  10cm  diameter 
plugs  extracted  for  measurements  and  the  effects  of  disturbing  the  soil. 
Results  of  these  measurements  are  shown  in  Table  1.  The  density  of  the 


TABLE  1 


Measured  Values  of  Flow  Resistance 

Site 
Sardis 

UM  Test  Field 

nonporous  lower  layer  material  at  the  UM  test  field  was  determined  from  a 

sample  of  clay  extracted  from  approximately  0.5  meters  below  the  surface. 

At  Sardis  beach,  the  density  of  the  clay  underlying  the  sand  was  assumed 
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to  be  the  same  as  at  the  UM  test  field,  which  was  found  to  be  2.65g/cm  . 

The  shear  and  compressional  wave  speeds  and  the  depth  of  the  porous 
material  were  determined  from  a  ray  propagation  model ^  for  a  two  layer 
system,  see  Figure  3.6.  An  impulse  is  created  at  x^  by  striking  a  metal 
block  with  a  hammer.  A  seismic  sensor  is  located  at,  say  x^  and  is  sensi¬ 
tive  to  both  shear  and  compressional  waves,  and  it  is  assumed  that  there 
are  three  possible  paths  to  the  sensor.  Path  one  is  along  the  surface  at  a 
wave  speed  of  v^  and  is  referred  to  as  the  direct  path.  The  second  path  or 
reflected  path  is  down  towards  the  lower  boundary  and  reflected  in  accor¬ 
dance  with  Snell's  law.  The  third,  the  refracted  path,  is  down  towards  the 
lower  boundary,  along  the  surface  of  the  lower  medium  at  a  speed  v^  then 
back  towards  the  upper  boundary.  When  the  seismic  sensor  is  close  to  the 
source  compared  to  the  layer  depth,  the  first  arrival  will  be  the  direct 
pulse.  As  the  seismic  sensor  is  moved  away  from  the  source  and  the  separa¬ 
tion  distance  becomes  greater  than  twice  the  layer  depth,  the  first  arrival 
will  be  the  refracted  pulse.  A  plot  of  arrival  times  versus  source  detec¬ 
tor  separation  distance  for  the  UM  test  field  is  plotted  in  Figure  3.7  and 
the  wave  velocities  are  the  slopes  of  the  lines  drawn.  From  the  geometry 
of  the  model,  the  layer  depth  can  be  expressed  as 
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where  x  is  the  location  of  the  intersection  of  the  lines  on  the  time- 
c 

position  graph.  Physically,  xc  is  the  distance  from  the  source  to  the 
point  at  which  the  refracted  pulse  arrives  at  the  same  time  as  the  direct 
pulse.  The  results  of  the  wave  speeds  and  the  layer  depths  for  both  sites 
are  recorded  in  Table  2. 
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TABLE  2 


Seismic 

:  Wave  Velocities 

(cm/s) 

t 

Layer  depth  (cm) 

$ 

Sardis  Beach 

s 

First  Layer 

Second  Layer 

Jj 

hi 

Compressional  Shear 

Layer  Depth 

Compressional  Shear 

S 

1.6xl04  1.5xl04 

26 

4.6xl04  2.6xl04 

UM  Test  Field 

First  Layer 

Second  Layer 

Compressional  Shear 

Layer  Depth 

Compressional 

Shear 

1.6xl04  . 95xl04 

30 

3.7xl04 

2.2x10 

Sources  and  Sensors 

The  sound  source  used  to  broadcast  acoustic  signals  was  constructed  at 
the  University  of  Mississippi  Physics  Department  and  a  diagram  of  the 
speaker  system  is  shown  in  Figure  3.8.  The  outer  most  material  is  1/2  inch 
particle  board  which  was  nailed  and  glued  to  an  internal  framework 
constructed  from  two  inch  by  two  inch  lumber.  The  speaker  system  consists 
of  four  independent  vented  boxes  of  which  each  has  an  internal  volume  of 
725  liters  and  are  each  tuned  to  20  hertz  with  a  four  inch  diameter  by  six 


32 


vH'  O  C  *w*‘o‘0* 


v  » bvv,  vs  y. 


inch  depth  cylindrical  port.  On  the  front  face  of  each  box  is  a  modified 
Peavy  18  inch  Black  Widow  driver.  The  speaker  system  was  driven  by  two  ADS 
Power  Plate  100  watt  amplifiers  and  was  mounted  on  a  truck  so  that  it  would 
be  made  mobile  for  field  use. 

A  Hall  Engineering  acoustic  signal  generator  was  used  to  record  two 


minute  segments  of  pink  noise  at  each  one  of  the  third  octave  bands  from  25 
Hz  to  500  HZ.  In  the  field  experiments,  this  tape  was  broadcast  using  a 
Nakamichi  550  tape  player.  A  block  diagram  of  the  acoustic  broadcasting 
network  is  shown  in  Figure  3.9. 

The  seismic  sensor  used  was  a  commercially  available  Mark  Products  L- 

4-3-D  triaxial  geophone.  This  geophone  contains  three  coils  which  are 

sensitive  to  motion  in  the  transverse,  radial  and  vertical  directions.  A 

calibration  factor  supplied  by  the  manufacturer  enables  the  determination 

of  the  seismic  velocity  from  the  signal  voltage.  The  calibration  factor 

-4 

was  1.0  mV  =  3.57  x  10  cm/sec.  Using  an  auger,  a  cylindrical  hole  was  dug 
slightly  larger  than  the  diameter  of  the  geophone.  The  geophone  was  then 
placed  in  the  hole,  leveled,  and  then  covered  with  the  soil  or  sand  that 
had  been  removed.  In  the  transfer  function  experiment  the  geophones  were 
buried  to  depths  such  that  the  top  surface  of  the  geophone  was  either  5  or 
25  cm  below  the  ground. 

The  probe  microphone  was  of  our  own  design  and  a  cut-away  view  of  the 
instrument  is  shown  in  Figure  3.1.  At  the  UM  test  site  the  probe  depth  was 


8cm  and  at  Sardis  beach  the  depth  was  12cra.  A  reference  AkG  microphone, 
identical  to  the  one  in  the  probe,  was  positioned  approximately  5  cm  above 
the  surface  of  the  ground. 

Figure  3.10  shows  the  relative  positions  of  the  sensors  and  sound 
source.  Four  different  angles  of  incidence,  in  increments  of  5°,  were  used 
during  the  experiment.  At  the  largest  angle,  20°,  the  separation  distance 
between  the  source  and  sensors  was  5.5  m  while  at  the  smallest  incident 
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angle  it  was  23  m.  The  angle  of  incidence  was  changed  by  repositioning  the 
truck  upon  which  the  speakers  were  mounted. 

Each  of  the  electronic  signals  from  the  sensors  were  amplfied  with  a 
Tektronix  AM  502  differential  amplifier  and  then  recorded  on  a  Teac  R-81 
seven  channel  FM  cassette  recorder.  A  block  diagram  of  the  data  recording 
system  is  shown  in  Figure  3.11.  The  amplifier  gains  were  adjusted  such 
that  the  most  intense  sound  level  broadcast  did  not  overload  the  ampli¬ 
fiers.  The  gains  of  each  of  the  coils  of  the  geophone  and  the  microphones 
along  with  the  recorder  channel  assignments  are  in  Table  3.  The  output 
voltage  of 


TABLE  3 
UM  Test  Field 
SENSOR 

Vertical  Triax 
Longitudinal  Triax 
Transverse  Triax 
Above  Ground  Microphone 
Probe  Microphone 
Voice  Memo  &  DC  Trigger 


Sardis 

SENSOR 

Vertical  Traix 
Transverse  Triax 
Longitudinal  Triax 
Above  Ground  Microphone 
Probe  Microphone 
Voice  Memo  &  DC  Trigger 


each  amplifier  was  five  volts  and  that  was  the  maximum  input  voltage  of  the 
tape  recorder,  thus  allowing  the  maximum  dynamic  range  of  the  recorder,  40 
dB,  to  be  used.  After  the  gains  were  set,  a  General  Radio  model  1986  sound 


level  calibrator  was  used  to  record  a  94  dB  calibration  tone  on  the  refer¬ 
ence  microphone.  The  output  voltage  of  this  reference  microphone  was 
measured  and  then  temporarily  connected  to  each  input  of  the  geophone 


Sound  Source 


Figure  3. 10.  Diagram  of  source  and  sensor  geometry  indicating  source 
height  h  ,  reference  microphone  position  h  ,  probe  depth  d  and  geophone 
depth  d  .s  m  p 

g 


Figure  3.11.  Data  recording  system  used  in  the  measurements  of  the  trans¬ 
fer  functions. 


signal  amplifiers,  thus  recording  a  calibration  voltage  for  the  seismic 
sensors.  This  procedure  for  calibration  tones  was  repeated  after  the 
experiment  was  complete.  Once  the  initial  calibration  tones  were  recorded, 
the  test  tape  containing  one-third  octave  bands  of  pink  noise  was  broad¬ 
cast.  The  geophone,  probe  and  reference  microphone  responses  were  recorded 
on  the  assigned  tape  recorder  channels.  A  five  second  1.5  volt  d.c.  pulse 
was  also  recorded  on  its  assigned  channel  (see  Table  3)  at  the  beginning  of 
the  broadcast  of  each  one  1/3  octave  band.  This  pulse  was  used  as  a 
trigger  when  the  data  were  digitized.  After  the  last  1/3  octave  band  was 
broadcast,  the  source  tape  was  rewound,  the  truck  repositioned  and  the 
broadcast  repeated.  This  proceedure  was  repeated  at  each  of  the  remaining 
incident  angles. 
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4.  DATA  ANALYSIS 


The  primary  quantity  of  interest  in  this  study  was  the  transfer  func¬ 
tion  defined  by 

Signal  Below  the  Surface  (Complex) 

T(f )  -  - . 

Sound  Pressure  Above  the  Surface  (Complex) 

The  transfer  function,  T(f),  is  a  complex  quantity  which  can  be  expressed 
in  terms  of  a  magnitude  and  phase.  The  data  analysis  procedure  was 

designed  to  extract,  from  the  FM  tape  records  made  in  the  field,  this 
magnitude  and  phase. 

Since  we  are  concerned  with  the  absolute  magnitude  of  the  transfer 
function,  the  first  step  in  the  analysis  was  to  determine  calibration 
levels  of  the  two  channels  to  be  compared.  Calibration  tones  were  gene¬ 
rated  for  each  data  channel  used  (except  the  probe  microphone  see 
experimental  procedures).  With  the  geophone  and  calibration  levels  used,  a 
1  bar  acoustic  sound  pressure  level  and  a  geophone  velocity  of 

1.754  x  10  ^cm/sec  produced  equal  voltage  levels  at  the  input  to  the  FM 
recorder. 

The  field  records  were  digitized  using  a  PDP  11/23  contolled  A/D 

converter  digitizing  at  1000  Hz.  No  anti-aliasing  filter  was  used.  We 

assumed  that  the  electronic  and  geophone  frequency  roll  was  sufficient  to 
prevent  aliasing. 

The  program  used  to  analyze  the  data,  called  Tranff.Bas,  is  attached 
as  Appendix  A. 7. 

Representative  results  of  the  data  analysis  are  given  in  Figures  4.1, 
4.2,  4.3,  and  4.4.  The  total  amount  of  data  precluded  a  complete  presen¬ 

tation  here  but  all  data  will  be  provided  on  request.  A  larger  sample  is 
included  in  Appendix  B.  Referring  to  Figure  4.2  note  that  the  direction 
of  phase  change  with  frequency  depends  upon  the  sign  convention  chosen  on 
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Figure  4 
Cion  at 


•  1.  The  experimental  magnitude  of  the  vertical  seismic  transfe 
at  a  depth  of  5cm  and  20°  angle  of  incidence  at  Sardis. 


2.  The  experimental  phase  of  the  vertical  seismic  transfer  func 
depth  of  5cm  and  10°  angle  of  incidence  for  the  UM  test  field. 


exp(iuit).  Unfortunately,  all  data  analysis  was  done  independently  of  the 
theory  development  so  the  signs  are  reversed.  In  the  data  analysis  pro¬ 
gram,  the  phase  becomes  more  negative  as  the  frequency  increases  (fixed 
path  length  difference,  decreasing  wavelength);  in  the  theory  program,  this 
would  result  in  a  more  positive  phase.  Also  note,  in  figure  4.3  that  the 
probe  had  no  absolute  calibration  so  we  only  measured  the  change  in  probe 
response  as  a  function  of  frequency.  For  convenience  in  presentation,  it 
was  assumed  that  the  probe  microphone  and  the  AKG  microphone  in  the  probe 
had  the  same  sensitivity  but  this  could  be  off  by  as  much  as  a  factor  of 
two.  Finally,  note  that  all  phases  are  limited  between  +90°  due  to  the 
ambiguity  in  the  arctan  routine.  This  accounts  for  occasional  discontin¬ 
uous  jumps  from  -90°  to  +90°.  In  fact,  the  phase  is  smoothly  varying  at 
this  point  through  an  angle  greater  than  -90°. 

Probe  microphone  data  collected  in  the  UM  test  field  suffered  from 
poor  S/N  due  to  the  rapid  attenuation  of  the  sound  field  with  depth  at  that 
site.  The  poor  S/N  ratio  makes  the  phase  measurements  especially  ragged 
through  the  trend  is  clear. 


5.  COMPARISON  OF  EXPERIMENTAL  DATA  TO  THEORETICAL  CALCULATION 


This  section  provides  a  physical  interpretation  of  selected  experimen¬ 
tal  data  and  a  comparison  to  theoretical  predictions.  The  effects  of  the 
parameters  which  characterize  the  porous  layer,  viscosity,  shape  factor 
ratio,  tortuosity,  flow  resistivity,  porosity,  angle  of  incidence,  and 
layer  depth,  on  the  acoustic  and  seismic  transfer  functions  are  considered. 
Also,  using  measured  values  of  flow  resistivity,  an  acoustic  surface  impe¬ 
dance  is  calculated  using  a  standared  empirical  form  and  compared  to  that 
which  the  Biot-Stoll  Model  (developed  here)  predicts  and  experiment. 


Experimental  Observations 


At  either  depth  at  the  UM  test  field,  the  phase  of  the  seismic 
transfer  function  shows  only  a  slight  change  when  the  incident  angle  is 
changed  from  10°  to  20°  for  frequencies  above  100  Hz  (See  Figure  5.1  and 
5.2).  Similar  observations  can  be  made  for  the  phase  measurements  at 
Sardis  beach,  but  with  less  assurance  due  to  the  more  rapid  oscillation  of 
the  phase  frequency  (Figure  5.3  and  5.4).  The  phase  of  the  seismic 
transfer  function  at  the  UM  test  field  shows  a  -|-Tr  phase  change  over  300  Hz 
frequency  range  (Figure  5.1);  at  Sardis,  for  the  same  depth,  there  is  a  3tt 
to  4ir  phase  change  for  the  same  frequency  span  (Figure  5.3).  At  a 
frequency  of  100  Hz,  there  is  approximately  a  30°  phase  difference  between 
the  5  cm  and  25  cm  depths  for  the  seismic  transfer  functions  at  the  UM  test 
field,  (Figure  5.1  and  5.4)  (angle  of  incidence  20°)  which  would  correspond 
to  a  seismic  phase  velocity  of  2.4  x  10  cm/s.  Table  2  shows  values  of 
measured  wave  speeds  near  this  value. 

The  magnitude  of  the  seismic  transfer  function  is  insensitive  to  angle 
of  incidence  at  both  depths  at  each  experimental  site.  However,  the  magni¬ 
tude  of  the  acoustic  transfer  function  at  Sardis  at  the  12  cm  depth  shows 
an  increase  with  decreasing  angle  of  incidence.  At  100  Hz,  there  is  a  40% 
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Figure  5.1.  The  phase  of  the  experimental  seismic  transfer  function  at 
5cm  depth  at  the  UM  test  field  for  10°  ( - )  and  20°  ( - ). 


Figure  5.2.  The  phase  of  the  experimental  seismic  transfer  function  at 
25cm  depth  at  the  UM  test  field  for  10°  ( - )  and  20°  ( - ). 
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Figure  5.3.  The  phase  of  experimental  seismic  transfer  function  at  a  5cm 
depth  at  Sardis  beach  for  10°  ( - )  and  20°  (— ). 
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Figure  5.4.  The  phase  of  experimental  transfer  function  at  a  25cm  depth  at 
Sardis  beach  for  10°  ( — — )  and  20°  ( - ). 


p 


«•- 


.*• 


.  i 

.v 

ci 


factor  ratio,  0.4iJ=<1.0,  or  tortuosity,  1.0<q<4.0,  the  acoustic  or  seismic 
transfer  functions  are  changed  by  less  than  2%.  The  magnitude  of  the 
transfer  functions  are  also  insensitive  to  changes  in  both  dynamic  and 
kinematic  viscosities  of  air  when  changed  by  factors  up  to  five.  The  flow 
resistance  and  porosity  do  have  a  significant  effect  on  the  transfer 
functions. 

The  effect  of  increasing  the  porosity  from  .3  to  .5  on  the  magnitude 
of  the  seismic  transfer  function  is  shown  in  Figure  5.6.  Increasing 
porosity  results  in  an  increase  in  the  seismic  transfer  function;  the 
magnitude  of  the  acoustic  transfer  function  decreases  for  the  same 
increases  in  porosity  as  indicated  in  Figure  5.7.  The  increase  in  the 
seismic  transfer  function  results  from  a  decrease  in  the  mass  density  of 
the  frame.  For  the  same  amount  of  incident  acoustic  energy  when  more  is 
coupled  into  the  frame,  less  must  be  available  as  transmitted  acoustic 
energy,  consequently  the  acoustic  transfer  function  decreases  with  the 
increase  in  porosity.  In  the  real  world,  the  porosity  cannot  change  with¬ 
out  changing  the  flow  resistivity.  The  argument  here  is  offered  only  to 
aid  physical  interpretation. 

The  effect  of  changing  the  flow  resistivity  from  60  rayls/cm  to  420 
rayls/cm  to  3000  rayls/cm  on  the  magnitude  of  the  seismic  transfer  func- 


tions  is 

shown 

in  Figure  5.8. 

Figure  5.9  shows 

the 

magnitude 

of 

the 

acoustic 

transfer 

function  for 

flow 

resistivities 

of 

60  rayls/cm. 

150 

rayls/cm 

and  420 

rayls/cm. 

The 

magnitude  of 

the 

acoustic 

transfer 

function  decreases  as  the  flow  resistivity  is  increased.  The  magnitude  of 
the  seismic  transfer  function  decreases  with  increasing  flow  resistivity 
for  frequencies  above  130  Hz  and  increases  at  the  lower  frequencies.  If  the 
flow  resistivity  were  to  go  to  infinity,  the  magnitude  of  the  acoustic 
transfer  function  would  tend  to  zero  because  acoustic  energy  could  not  be 
transmitted  across  the  boundary  as  the  pore  fluid  could  not  move.  The 
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Figure  5.6.  Magnitude  of  the  calculated  seismic  transfer  function  for 
porosities  of  .3  (— — )  and  .5  ( — -) ,  5  cm  below  the  surface  and  5  angle  of 
incidence.  All  other  parameters  are  set  to  values  measured  at  Sardis. 
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Figure  5.7.  Magnitude  of  the  calculated  acoustic  transfer  function  for 

porosities  of  .3  ( - )  and  .5  ( - ),  5  cm  below  the  surface  and  5°  angle  of 

incidence.  All  other  parameters  are  set  to  the  values  measured  at  Sardis. 
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Figure  5.8.  Magnitude  of  the  calculated  seismic  transfer  function  for  flow 
resistivities  of  60  rayls/cm  (— — ),  420  rayles/cm  (-•-)  and  3000  rayls/cm 

( - )  at  23cm  and  20  .  The  porosity  is  .4  and  other  parameters  are  set  tr 

the  values  measured  at  Sardis. 
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Figure  5.9.  Magnitude  of  the  calculated  acoustic  transfer  function  at  5cm 

for  flow  resistivities  of  60  rayls/cm  ( - ),  150  rayls/cm  £ - )  and 

420  rayls/cm  (-•-),  a  porosity  of  .4  and  angle  of  incidence  of  5°.  Other 

parameters  are  set  to  the  values  measured  at  Sardis. 
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seismic  transfer  function  would  approach  zero  since  no  acoustic  energy 


could  be  coupled  to  motion  of  the  frame  through  drag  forces  at  the  pore 


walls.  If  the  flow  resistivity  went  to  zero,  all  the  incident  acoustic 


energy  would  pass  through  the  upper  boundary  without  attenuation  but  there 


would  again  be  no  seismic  motion  since  there  would  be  no  frame. 


By  changing  the  layer  depth,  the  effect  of  the  measured  wave  speeds 


and  consequently  the  elastic  moduli  on  the  transfer  functions  can  be 


observed.  Figure  5.10  shows  the  seismic  transfer  function  for  layer  depths 
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Figure  5.10.  Magnitude  of  the  calculated  seismic  transfer  function  at  23cm 

for  layer  depths  of  28  cm  (  — )  and  24  cm  ( - ).  Angle  of  incidence, 

porosity  and  flow  resistance  are  5  ,  .4,  350  rayls/cm.  Other  parameters 
are  set  to  those  measured  at  Sardis. 


of  28  and  24  cm.  The  acoustic  transfer  function  for  the  same  layer 


depths,  not  shown,  is  changed  by  less  than  12.  The  acoustic  transfer 


function  is  insensitive  to  these  changes  but  the  seismic  transfer  function 


shows  frequency  shifts  in  the  minimum  and  maximum  indicating  that  some  of 


the  seismic  energy  is  reflected  from  the  lower  boundary. 


,VvV.- 


During  the  transfer  function  experiments,  the  angle  of  grazing 
incidence  was  varied  from  5°  to  20°.  Figures  5.11  and  5.12  indicate  the 
calculated  effect  of  angle  of  incidence  on  the  magnitude  and  phase  of  the 
seismic  transfer  function.  The  magnitude  and  phase  of  the  acoustic 
transfer  function  is  insensitive  to  the  angle  of  incidence  to  less  than  1 
percent  while  the  magnitude  of  the  seismic  transfer  function  decreases  by 
as  much  as  25%  at  the  maximum  and  the  phase  curve  is  shifted  vertically  by 
4  to  5  . 

Acoustic  Transfer  Function 

Next,  the  experimental  acoustic  transfer  function  will  be  compared  to 
the  predicted  transfer  function  using  the  measured  values  of  flow 
resistivity  (See  Table  1)  with  the  porosity  set  to  0.4  for  both  experimen¬ 
tal  sites.  Figures  5.13  and  5.14  show  the  magnitude  and  phase  of  both 
the  experimental  and  predicted  acoustic  transfer  functions  for  the  Sardis 
beach  area  at  an  incident  angle  of  20°  and  12  cm  below  the  surface.  The 
theoretical  phase  plotted  here  is  the  negative  of  what  was  predicted  for 
reasons  previously  explained.  Comparisons  of  transfer  functions  for  the  UM 
test  field  are  not  shown.  The  agreement  is  no  better  or  worse  than  at 
Sardis. 

Even  when  the  layer  depth  is  theoretically  changed  from  26  cm  to  6  cm 
the  magnitude  of  the  transfer  function  5  cm  below  the  surface  does  not  show 
interference  with  waves  reflected  from  the  lower  boundary  because  the 
acoustic  wave  in  the  pores  is  rapidly  attenuated;  it  is  not  reflected  from 
the  lower  boundary.  By  the  time  that  the  acoustic  wave  has  propagated  4  cm 
through  the  pores  almost  all  the  acoustic  energy  has  been  transferred  to 
motion  of  the  frame  at  the  pore  walls.  The  effect  of  the  porous  layer 
acting  as  an  impedance  matching  device  coupling  acoustic  energy  into  the 
non-porous  clay  will  be  demonstrated  in  the  next  section  of  this  report. 
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re  5.11.  The  calculated  effect  of  angle  of  incidence  on  the  magnitude 
the  seismic  transfer  function  at  5°  (— )  and  20  ( - )  at  a  depth  of 


Figure  5.11. 
of  the  seis 
23cm.  The 
Sardis. 


The  porosity  is  0.4  and  other  parameters  are  the  values  measured  at 
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Figure  5.12.  The  calculated  effect  of  angle  of  incidence  on  the  phase  of 

the  seismic  transfer  function  at  5  ( - )  and  20  ( - )  at  a  depth  of  23cm. 

The  porosity  is  0.4  and  other  parameters  are  the  values  measured  at  Sardis. 


Seismic  Transfer  Function 


In  the  transfer  function  experiment,  the  seismic  sensor  was  buried  so 
that  the  top  of  the  case  was  at  depths  of  5  and  25  cm.  However,  because 
of  the  physical  size  of  the  sensor  (cylindrical,  diameter  =  20  cm,  length  = 
18  cm),  the  transfer  functions  are  calculated  at  23  cm  and  43  cm  below  the 
surface.  At  a  depth  of  43  cm,  the  lower  clay  layer  is  reached  so  we 
measure  the  seismic  transfer  function  in  both  the  porous  soil  and  the  non- 
porous  clay.  Figure  5.15,  5.16,  5.17  and  5.18  show  the  comparison  of  the 
magnitude  and  phase  of  the  experimental  and  predicted  seismic  transfer 
function  for  both  depths  at  Sardis.  The  experimental  magnitude  shows  a 
strong  minimum  at  160  Hz  and  other  fine  structure  which  is  not  predicted  by 
our  single  layer  model.  At  Sardis  the  moisture  content  of  the  sand 
increased  rapidly  below  5  cm  and  became  saturated  at  the  depth  of  the  clay 
layer.  To  predict  the  detailed  structure  of  the  seismic  transfer  functions 
in  such  a  media  it  may  be  necessary  to  incorporate  multiple  layers  in  the 
physical  model.  In  Figure  17,  the  UM  test  field  measured  seismic  wave 
speeds  for  the  clay  were  also  used  in  the  calculations  along  with  a  layer 
depth  of  35  cm.  Comparisons  of  the  UM  test  field  data  to  theory  are  not 
shown,  but  again  the  agreement  is  as  good  as  at  Sardis. 

From  the  measured  seismic  wave  speeds  in  the  clay  below  the  porous 
medium,  the  reflection  coefficient  for  an  air-clay  boundary  would  be 
approximately  .999  indicating  almost  no  transmission  into  the  surface. 
However,  the  measured  magnitude  of  the  seismic  transfer  function  for  the 
lower  medium  indicates  as  much  seismic  energy  is  transmitted  into  the  clay 
as  is  coupled  into  the  porous  medium  from  the  air.  Theoretically,  this  is 
also  seen.  Transmission  coefficients  of  .75  are  computed  at  200  Hz 


indicating  that  the  porous  layer  acts  as  an  impedance  matching  device  for 
coupling  acoustic  energy  into  the  non-porous  clay  below  the  porous  medium. 
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Figure  5.15.  The  magnitude  of  calculated  ( - )  and  experimental  ( 

seismic  transfer  function  at  Sardis  at  20  and  23  cm. 
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Figure  5.16.  The  phase  of  calculated  ( - )  and  experimental  ( - )  seismic 

transfer  function  at  Sardis  at  20°  and  23  cm. 
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Figure  5.17.  The  magnitude  of  the  calculated  (— •— )  and  experimental  (— —  ) 
seismic  transfer  function  at  Sardis  at  20°  and  43  cm.  The  solid  line  is 
calculated  with  the  UM  test  field  clay  velocities  and  a  layer  depth  of  35cm 
substituted  for  the  Sardis  values* 


Figure  5.18.  The  phase  of  the  calculated  ( - )  and  experimental  ( - ) 

seismic  transfer  function  at  Sardis  at  20°  and  43  cm. 


The  very  fact  that  the  theoretical  predictions  come  within  a  factor  of  2  of 
the  experimental  data  in  all  cases,  compared  to  a  factor  of  1000  for  a 
homogeneous  layer  theory,  is  remarkable  considering  that  there  are  no 
empirical  adjustments  or  parameters  allowed. 


Surface  Impedance 


The  acoustic  surface  impedance  at  Sardis  or  the  UM  test  field  soil  was 
not  experimentally  determined.  However,  Bolen  and  Bass  have  published 
normalized  surface  impedances  for  several  ground  covers  which  include  a 
sand  and  two  different  soils.  They  compared  experimentally  measured  values 
of  surface  impedance  to  that  which  is  predicted  from  a  single  parameter 
model  based  on  measured  flow  resistivity  (see  Reference  17).  In  each  of 
the  ground  covers  investigated,  they  found  that  the  impedance  necessary  to 
reproduce  the  sound  pressure  level  measured  with  above  ground  microphones 
required  that  a  they  use  a  flow  resistivity  in  the  single  parameter  model 
approximately  one-half  the  measured  flow  resistivity  (see  Figure  5.19, 
5.20).  This  was  referred  to  as  the  best  fit  to  measured  values.  The 
theoretical  normalized  surface  impedance  predicted  by  the  Biot-Stoll  Model 
developed  here  is  plotted  on  the  Bolen  and  Bass  impedance  curves  for  their 
measured  values  of  flow  resistivity.  Note  that  no  empirical  adjustment  of 
the  flow  resistance  is  required.  In  fact,  there  are  no  adjustable 
parameters  allowed. 
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Figure  5.19.  The  real  and  imaginary  impedance  values  for  a  Florida  Sand 
determined  by  Bolen  and  Bass;  experimental  data  (x),  empirical  Chessell 

model  ( - ),  the  best  fit  ( - ).  The  dash  dot  line  is  calculated  using  the 

Biot-Stoll  Model  with  the  measured  values  of  flow  resistance, 
o  ■  60  rayls/cm. 
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Figure  5.20.  The  real  and  imaginary  impedance  values  for  a  Softball  field 

determined  by  Bolen  and  Bass;  experimental  data  (x),  empirical  Chessell 

model  ( - ),  the  best  fit  ( - ).  The  dash  dot  line  is  calculated  using  the 

Biot-Stoll  Model  with  the  measured  value  of  flow  resistance, 

o  ■  400  rayls/cm. 
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6.  SUMMARY  AND  CONCLUSION 


Using  the  Biot-Stoll  Model  of  a  porous  elastic  medium  the 
compressional  wave  propagation  constants  have  been  numerically  determined 
for  a  porous  sand  and  a  grass  covered  field.  The  earth's  surface  was 
modeled  as  a  finite  layer  of  porous  material  between  two  semi-infinite 
media,  air  and  non-porous  clay.  The  complex  displacement  amplitudes  of 
compressional  and  shear  waves  in  each  medium  were  found  using  numerical 
techniques  and  a  complex  transfer  function  was  defined  in  terms  of  these 
amplitudes.  Compressional  waves  are  emphasized  here  to  clarify  the  physi¬ 
cal  properties.  The  surface  impedance  was  also  determined  from  the  wave 
amplitudes.  An  outdoor  experiment  was  performed  in  order  to  measure  the 
transfer  function  so  that  theory  could  be  compared  to  experiment. 

A  sensitivity  analysis  of  the  calculated  transfer  function  revealed 
that  the  only  physical  quantities  which  significantly  affect  the  computed 
transfer  funtions  are  the  flow  resistivity  of  the  porous  surface,  the  depth 
of  the  porous  layer,  and  the  compressional  and  shear  wave  velocities  of  the 
underlying  clay.  The  magnitudes  of  the  calculated  acoustic/seisraic 
transfer  function  using  independently  measured  values  of  the  physical 
properties  typically  agrees  with  the  measured  values  of  the  transfer  func¬ 
tion  within  factors  of  two.  The  theory  does  not  predict  observed  fine 
structure  in  the  acoutic/seismic  transfer  function.  It  may  be  necessary  to 
include  multiple  layers  and  flow  resistivity  gradients  in  the  porous 
surface  if  the  fine  structure  is  required. 

The  Biot-Stoll  Model  allows  for  an  acoustic  wave  in  the  pores  of  the 
porous  surface.  In  a  separate  experiment,  attenuation  coefficients  and 
phase  velocities  of  this  acoustic  wave  were  measured  in  sand  using  a  probe 
microphone.  The  phase  velocity  of  the  acoustic  wave  in  the  pores  was  only 
tens  of  m/s  and  the  attenuation  was  one  to  two  dB/cm.  Predicted  attenua¬ 
tion  and  wave  velocity  agree  well  with  experiment. 


The  mechanism  for  the  coupling  of  acoustic  energy  into  the  ground  as 
seismic  energy  is  understood  in  the  light  of  these  measurements.  As  the 
airborne  acoustic  wave  impinges  upon  the  porous  surface,  it  forces  motion 
of  air  in  the  surface  pores.  The  oscillating  air  in  the  pores  is  subjected 
to  viscous  drag  at  the  pore  walls.  As  a  result  of  this  drag,  energy  is 
taken  out  pore  fluid  motion  and  appears  as  motion  of  the  pore  walls  - 
motion  of  the  soil  matrix.  This  motion  is  detected  with  a  seismic  sensor  - 
a  geophone.  An  alternative  way  of  viewing  this  coupling  process  is  to 
consider  the  porous  layer  as  an  impedance  matching  network.  Without  the 
pores,  the  acoustic  resistance  (pc)  in  air  is  much  less  than  in  the  earth. 
The  pore  wave,  however,  is  very  slow  bringing  the  effective  pc  product  for 
the  earth's  surface  more  near  that  of  air. 

This  physical  model  also  explains  a  result  which  has  mystified  outdoor 
acousticians  for  a  number  of  years.  Specifically,  when  one  measures  the 
impedance  of  the  earth's  surface  it  is  found,  in  almost  all  cases,  to  be 
locally  reacting.  This  means  that  energy  is  refracted  so  that  it  is  almost 
normal  to  the  surface.  The  impedance  then  becomes  essentially  independent 
of  the  angle  of  incidence  and  regions  of  the  surface  very  far  away  from  the 
point  where  the  acoustic  wave  strikes  do  not  contribute  to  the  refracted  or 
reflected  acoustic  fields.  This  is  difficult  to  understand  if  the  wave 
speed  in  the  earth  is  taken  to  be  hundreds  of  meters  per  second  but  follows 
readily  from  Snell's  Law  if  the  wave  speed  in  the  surface  (in  the  pores)  is 
only  tens  of  meters  per  second.  The  empirical  result  that  the  earth's 
surface  is  almost  always  locally  reacting  can,  then,  be  understood  as  a 
natural  consequence  of  the  porous  upper  layer  of  the  earth. 

In  conclusion,  the  Biot-Stoll  Model  for  a  porous  elastic  medium  as 
presented  here  provides  a  good  physical  description  of  the  coupling  of 
airborne  acoustic  energy  into  the  ground  where  it  appears  as  seismic  motion 
and  an  acoustic  wave  in  the  pores.  This  model  predicts  the  measured 


attenuation  of  the  acoustic  wave  in  the  pores,  surface  impedance,  and 
acoustic  and  seismic  transfer  functions  and  their  dependence  upon  angle  of 
incidence  and  frequency  within  factors  of  two  of  measured  values  even  when 
one  assumes  only  one  porous  homogeneous  layer.  Earlier  calculation  based 
upon  semi-infinite  homogeneous  fluids  in  contact  were  in  error  by  a  factor 
of  1000  and  these  calculations  predicted  a  strong  angle  of  incidence  not 
observed  experimentally,  and  predicted  incorrect  frequency  dependence. 

The  theoretical  results  presented  here  need  three  refinements  in  order 
to  make  this  procedure  a  useful  engineering  tool: 

1.  Only  plane  incident  waves  are  included  when,  in  fact,  most  sources 

18 

are  more  nearly  spherical.  Dr.  Attenborough  has  developed  the  mathema¬ 
tics  necessary  to  treat  the  spherical  case.  His  formalism  needs  to  be 
included  in  our  present  computer  programs  and  additional  experiments  need 
to  be  performed  for  angles  of  incidence  less  than  5°  where  Attenborough 
predicts  a  difference  between  plane  and  spherical  wave  results. 

2.  The  theoretical  expressions  need  to  be  simplified  by  searching  for 

realistic  approximations  so  that  one  can  treat  impulsive  sources  using 

Fourier  Transform  techniques.  This  refinement  will  require  extensive 

calculations  be  made  to  verify  that  the  approximate  relations  are  valid  for 

practical  surfaces.  This  refinement  will  allow  comparison  between  experi- 

19 

ment  and  the  large  data  base  developed  by  van  Hoof  and  co-workers  for 
impulsive  sources. 

3.  The  model  needs  to  be  extended  to  allow  for  variation  of  flow 
resistance  with  depth  and  multiple  porous  layers.  This  refinement  will  be 
necessary  before  we  can  hope  to  predict  the  fine  structure  in  the 
acoustic/seismic  transfer  functions.  Experimentally,  data  is  needed  for  a 
wider  variety  of  surfaces  as  well  as  data  with  incident  angles  less  than 
5°.  Specifically,  surfaces  with  very  low  and  very  high  flow  resistivities 
would  help  explore  the  range  of  validity  of  this  model. 
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The 

subroutine 

PROP. FORTRAN 

is  used  to  drive  the  routine  SOAIR. FORTRAN 

■ 

which  computes  the 

compressional  wave  propagation  constants  using  equation 

> 

1 

8.  The 

frequency, 

FR,  is  the 

input  parameter  and  the  propagation  con- 

? 

% 

0- 

stants , 

ZSLOW1,  ZFAST1  are  the  output.  The  routine  is  in  CGS  units.  SOAIR 

t 

computes 

the  viscosity  correction 

factor  (equation  6)  with  IARGT=0  and  the 

c 

k! 

J: 

complex 

compressibility  with  ICOMP-1  (see  Reference  12).  Parameters  not 

\ 

1 

defined 

in  SOAIR  are 

as  follows: 

* 

* 

NVISC 

Static  Shape  Factor  (n) 

1 

ft 

STATSF 

Dynamic  Shape  Factor  (s) 

v 

NTHER 

Same  as  NVISC 

W 

QSQ 

Tortuosity 

w 

r 

. 

RHO 

Bulk  Density 

> 

v 

CP 

Specific  Heat  of  Air 

< 

* 

COFUIS 

Same  as  ETA 

a 

J 

^  • 

THECON 

Thermal  Conductivity  of  Air 

* 

PRN 

Prandtl*s  Number 

•i \ 

SHFT 

Shape  Factor  Ratio 

• 

SHF 

Same  as  SHFT 

ft 

ft 

GAMA 

Ratio  of  Specific  Heats  of  Air 

T* 

COMP 

Compressibility  of  Air 

? 

tf 

HBAR 

See  equation  2a 

CBAR 

See  equation  2b 

m 

u 

V 

MBAR 

See  equation  2c 

k 

ft. 

D 

See  equation  2d 

t 

The 

subroutines 

CBESJ. FORTRAN  and  CBJ. FORTRAN  compute  the  Bessel 
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f 

.\ 

functions  required 

for  the  viscosity  correction  factor  and  complex 
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compressibility. 


DRIVE  SR  ROD TINE  FOR  PROPAGATION  PROGRAMS 
FOR  DEFINITION  OF'  CONSTANTS  SEE  APPROPRIATE 
SUBPROGRAMS.  SUBPROGRAM  SOAIR  COMPUTES  THE 
COMPLEX  PROPAGATION  CONSTANTS  FOR  THE  POROUS 
LAYER. 

SUBROUTINE  PROP(FR,ZSLCW1 ,ZFAST1 ) 

IMPLICIT  REAL  *8  (A-H.O-Z) 

REAL  *8  MU,  KB  BAR ,  K,  KR, NTHER,  NVISC,  N  PRIME 
COMPLEX  *16  ZSLOW1  ,ZFAST1  ,  J,  F 
COMMON  /BLK2 /  J,F,ETA,K,Q2 
COMMON  /BLK4/  PHI 

COMMON  /BLK/  DM, IARGT, ICOMP, MU, STATSF, NTHER, NVISC, NPRIME 
COMMON  /BLK3/  RHOF, GB, KBBAR,KR, RHOS, OMEGA 
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PR OOOOIO 
PR000020 
PR 000030 
PR 000040 
PR 000050 
PR000060 
PR 000070 
PR 000080 
PR 000090 
PR 000100 
PR 0001 10 
PR000120 
PR000130 
PR0001 40 


JiSV 


PR 000150 

ICOMP  =  0  PR 000160 
M0s.1495D0  PR 000170 
ETA=1 .81D-04  PR 0001 80 
K= ETA/PHI  PR000190 
CALL  SOAIB(FR,ZSLOW1 ,ZFAST1 )  PR000200 
RETURN  PR000210 
END  PR 000220 


StJB  ROD  TINE  SOAIR(FR, Z2L0W1  ,ZFAST1  ) 

INPUTS  FR 

OOTPOTS  ZS.0W1  ,ZFAST1 

TO  OBTAIN  PROPAGATION  CONSTANTS  FOR  DILITATIONAL  NAVES. 

15TH  NOV  82 

RECONSTRUCTED  VERSION,  25TH  APRIL  83. 

DATA  TAKEN  FROM  STOLL’S  PAPER  J ASA (1969)  V  47,  N  5, PI 440- 1447. 

SUBROUTINES  USED 

CBESJ 

CBJ 

SUBROUTINE  SOAXR(FR, ZSL0W1 , ZFAST1 ) 

IMPLICIT  REAL  «8  (A-H.O-Z) 

DIMENSION  AKBDB(6),FAT(6) 

DIMENSION  DEP(6),AT(6) 


REAL  *8  MU, KB BAR, K, KR, NTHER, NVISC, NPRIKE 

COMPLEX  «16  J, ROOTJ, CAPB, CAPC, SDISC, ARG, Z1 ,Z2 ,DISC 

COMPLEX  *16  XARG1 ,XJ0,XJ1 ,T1 tCOMP,KF 

COMPLEX  «16  D,EBAR,CBAR,MBARf CAPA 

COMPLEX  *16  XARG2,RJ0,RJ1,T2,M 

COMPLEX  *16  BJ0.BJ1 , T,F,BBYA,BBYA2 

COMPLEX  *16  ZKRB 

COMPLEX  *16  ZFAST1 ,AKB1 ,ZKRB1 ,ZSL0N1 

COMPLEX  »16  CPS, CPB, ZSLOW, ZFAST, SRAT1 , SRAT2 , FRAT1 , FRAT2 

COMPLEX  *16  XARG, AKB 

COMPLEX  «16  EP,EQ  ,ER,TCWS,TCWF,  RAFTOS 


COMMON  /BLK/  DM, IARGT, ICOMP, MU , STATSF, NTHER, NVISC, NPRIME 
COMMON  /BLK1/  EP,B3 ,ER 
COMMON  /BLK2/  J,  F,  ETA,  K,QS3 
COMMON  /BLK3/  RHOF,GB, KBBAR, KR, RHOS, OMEGA 
COMMON  /BLK7/  HBAR,CBAR,MBAR 


S0A00010 
SOA00020 
SOA00030 
SOA00040 
SOA00050 
SO  AO 0060 
SOA00070 
S0A00080 
SOA00090 
SOAOOIOO 
SOAOOIIO 
SOA00120 
SOA00130 
SOA00140 
SOA0015O 
SOA00160 
SOA00170 
SOA00180 
SOA00190 
SOA00200 
SOA0021 0 
SO AO 0220 
SOA00230 
SOA00240 
SOA00250 
SO AO 0260 
SOA00270 
SOA00280 
SOA00290 
SOA00300 
S0A00310 
S0A00320 
SOAO0330 
SOA00340 
S0A00350 
S0A00360 
SOA00370 
SO  AO 03  80 
SOA0039O 
SOA00400 
S0A00410 
SOA00420 
SOA0  0430 
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o  o  o  o  n  OOOOCIO  ooo  o  o  o  o  oooooooo  o  o 


PI  S  4 . D0*DATAN ( 1 . DO ) 

OMEGAs POROSITY( NO  UNITS). 

KRsBULK  MODULUS  OF  GRAINS  (DYNES/CM«2) 

MU = KINEMATIC  VISCOSITY  (CM”2/SEC) 

ETAs DYNAMIC  FLUID  VISCOSITY  (DYNE-SEC/CM»»2) 
RHOSs DENSITY  OF  SCLID  (GM/CM«3) 

RHOFsDENSITY  OF  FLUID (GM/CM»«3) 

RHO  =  ( 1 . DO-OMEGA) *RHOS  +  0MEGA*RH0F 


COMPUTE  PRANDTL  NUMBER. 

CP  s  1006. DO 

COFVIS  s  1.81D-05 
THECON  s  0.0257D0 

PRN  s  COFV IS  *CP/THEC0N 

SPRN  s  DSQRT(  PRN  ) 

Ks PERMEABILITY  OF  THE  POROUS  FRAME  (CM”2) 

SHFT  s  DSQRT(  STATS  FJ/NTHER 
SHF  s  DSQ  RT( STATSF) /NV ISC 


QSQ  3 

(OMEGA) ••C-NPRIME) 

GAMA 

s  1.4D0 

STWO 

3  DSQRT(2.D0) 

J 

3  (0.D0.1.D0) 

ROOTJ 

s  CDS3 RT(  J) 

ROOTJ  s  CMPLX(  REAL  (ROOTJ) ,-AIMAG(R00TJ) ) 
OM  s  2.D0*PI*FR 
0M2  s  OM*OM 
OM4  s  0M2*0M2 


COMPUTE  COMPLEX  COMPRESSIBILITY. 
XL  AMI  s  (8.DO«QS3*K)/(MU*OMEGA) 
XL AMI  s  SHFT*DSQ  RT( XL  AMI ) 

XL  AMI  s  XL AMI *DSQRT(0M) 

XARG1  s  ROOTJ »SPRN»XL AMI 
AXARG1 s  CDABS(XARGI) 

IF(  AXARG1  .GT.  7.0)  GO  TO  60 
CALL  CBESJ(XARG1 ,0,XJ0) 
CALL  CBESJ(XARG1 ,1 ,XJ1) 
60  IF(  AXARG1  .LE.  7.0)  GO  TO  70 
CALL  CBJ(XARG1 ,0,XJ0) 
CALL  CBJ(XARG1 ,1  ,XJ1 ) 

70  T1  s  XJ1/XJ0 


S0A00440 
S0A00450 
SOA0  0460 
S0A00470 
S0A00480 
SOA00490 
SOA00500 
SOA00510 
S0A00520 
SOA0  0  530 
S0A00540 
S0A00550 
SO  AO  0560 
SO  AO 0570 
SOA00580 
SOA00590 
SOA00600 
SOA00610 
SOA00620 
SOA00630 
SO AO 06 40 
SOA00650 
S0A00660 
SOAO0670 
SO AO 06 80 
SOA00690 
SOA00700 
SOA00710 
SOA00720 
SOA00730 
SOA00740 
S0A00750 
SOA00760 
S0A00770 
SO  AO  07  80 
S0A00790 
SOA00800 
SOA00810 
S0A00820 
S0A00830 
S0A00840 
S0A00850 
S0A00860 
SOA0O870 
S0A00880 
SOA00890 
SO AO 0900 
S0A0091 0 
S0A00920 
SOAO0930 
SOA00940 
SOA00950 
SO  AO  0960 
SOA00970 
S0A00980 
SOA00990 


r„  v.v.  ,v  V  >  v  vt- 


b  L 

L. 


a 

90 

•S' 

100 

110 

? 

V 

120 

.V 

,v, 

130 

c 

,% 

140 

C 

COMP  a  (1.D0+2.D0*(GAMA-1.D0)»T1/XARG1)/(GAMA«1.01325D06) 
KF  =  1. DO/COMP 
IF(  ICOMP  .NE.  1  }  GO  TO  80 
COMP  =  9.86923D-07 
KF  =  1. DO/ COMP 

CONTINUE 


D  a  KR*(1.D0  ♦  OMEGA*(KR/KF  -  1.00)) 

HBAH  =  (KR-KBBAR)»(KR-KBBAR)/(D-KBBAR) 

+  KBBAH  +  (4.00/3.D0)*GB 
CBAR  =  KR*(KH-KBBAR)/(D-KBBAR) 

MBAR  a  KR*KR/ ( D-KBBAR ) 

M  s  QSQ*(RH0F/0MEGA) 

CAPA  =  -HBAR*MBAR  +  CBAR*CBAR 

VISCOSITY  CORRECTION 

XARG  a  SHF*DSQRT(  (  8.DO*QS3*K)/(MU#OMEGA)  ) 

XX  =  XARG*DSQRT(OM) 

AEG  =  ROOTJ*XX 
0M2  a  OM*OM 
0M4  a  0M2*0M2 
IF(  IARGT  .NE.  1)  GO  TO  90 
ARG  a  ROOTJ*XX 
IF(  IARGT  . NE.  2)  GO  TO  100 

AP  *  SHF*DS3RT(  8.D0«QSQ«0MEGA) 

ARG  a  ROOTJ*DSQ RT( OM/ETA) *AP 
IF(  IARGT  .NE.  3)  GO  TO  110 

APHI  a  (DM/3. D0)»(0MEGA/(1. DO-OMEGA)) 
ARG  a  ROOTJ •APHI #DSQ  RT( OM/ETA ) 

IF(  IARGT  .NE.  4)  GO  TO  120 
ARG  a  ARG/OMEGA 

CONTINUE 

IF(  CDABS(ARG)  .GT.  7.0)  GO  TO  130 
CALL  CBESJ(ARGfO,BJO) 

CALL  CBESJ(ARG,1 ,BJ1) 

IF(  CDABS(ARG)  .LE.  7.0)  GO  TO  1 40 
CALL  CBJ(ARG«0 |BJ0) 

CALL  CBJ(ARG,1 ,BJ1) 

WRITE(6f •) *  BJO  a  *,BJ0,«BJ1  a  >,BJ1 
CONTINUE 
T  a  BJ1/BJ0 
WRITE(6,»)»  T  a  ' ,T 

F  a  (-0.25D0)»(ARG*T/(1 .DO  -  2.D0*T/ARG)) 
FREAL  a  DREAL(F) 

FIMAG  a  DIMAG (F) 

CAPB  a  0M2#(M#HBAR  -  2.D0*CBAR»RH0F  +  RHO*MBAR) 
♦J •OM*F*ETA*HB  AR/K 
CAPC  a  0M4  * ( RHOF*RH  OF  -  M*RHO) 

-  0M2*(RH0*J»0M»F»ETA/K) 

BBYA  a  -CAPB/ ( 2 . D0*CAPA) 

BBYA2  a  BBYA •BBYA 


SOAOIOOO 
SOAOIOIO 
SOA01020 
S0A01030 
SOA01040 
SOA01050 
S0A01060 
SOA01070 
SOA01080 
S0A01090 
SOAO 1100 
SOA01110 
SOAO 11 20 
SOAO 1130 
SOAO 11 40 
SOAO 11 50 
SOAO 1160 
SOAO 1170 
SOAO 11 80 
SOAO 11 90 
SOAO 1200 
SOAO 1210 
SOAO 1220 
SOAO 1230 
SOAO 1240 
SOAO 1250 
SOA0 1 26  0 
SOA0 1270 
SOAO 1280 
SOAO 1290 
SOAO 1300 
SOA01310 
SOAO 13 20 
SOAO 1330 
SOAO 13 40 
SOAO 1350 
SOAO 1360 
SOAO 1370 
SOAO 13 80 
SOAO 1390 
SOAO 1400 
SOA0 1410 
SOAO 1420 
SOA0 1430 
SOAO 1440 
SOAO 1450 
SOAO 1460 
SOAO 1470 
SOAO 14 80 
SOAO 1490 
SOAO 1500 
SOAO 1510 
SOAO 1520 
SOAO 1530 


V 

DISC  s  BBYA2  -  CAPC/CAPA 

SOA01540  ^ 

>  .. 

SDISCa  CDSQRT(DISC) 

SOA01550  v 

Z1  =  BBYA  ♦  SDISC 

S0A01560  ' 

c  *• 

Z2  s  BBYA  -  SDISC 

SOA01570  > 

u 

c  • ' 

Z1  a  CDSQRT(ZI) 

SOA01580  /. 

Z2  =  CDSQRT(Z2) 

SOA01590  | 

c 

SOA01600  % 

m 

ZSL0W1  a  Z1 

SO A0 1610  X 

ZFAST1  =  Z2 

S0A01620  £ 

Y  " 

c 

S0A01630 

L 

c 

SOA01640  j; 

1  'FI 

c 

S0A01650  | 

*s  .*  • 

RETURN 

SOA01660  ,V 

~  ’  * 

END 

SOA01670  < 

c 

SUBROUTINE  CBESJ(X,N,CBJ) 

CBB00010  Is 

y  - 

c 

INPUTS  X,N 

CBE00020  £ 

ISM 

c 

OUTPUTS  CBJ 

CBE00030  g. 

>  V* 

1/  h 

5  >j 

c 

ROUTINE  TO  COMPUTE  COMPLEX  BESSEL  FUNCTION  JN(X) 

CBE00040  < 

c 

WHERE  X  IS  A  COMPLEX  ARGUMENT 

CBE00050  j£ 

& 

3  v 

CBE00060  J 

c 

CBE00070  v; 

SUBROUTINE  CBESJ(X, N,  CBJ) 

CBE0  0  0  80  '£ 

IMPLICIT  REAL  »8  (A-H.O-Z) 

CBE00090  ' 

1  sM-ifcS? 

COMPLEX  *16  X,CBJ,FM1 ,FM, BMK,  ALPHA, CBPREV 

CBE00100 

P:  ,  :-: 

CBPREV  =  (0.0, 0.0) 

CBE00110  X 

A 

D  s  1.0E-4 

CBE00120 

P 

AX  a  CDABS(X) 

CBE00130  3 

ii 

NTEST  a  20.0  +  10.0*AX  -  AX»*2/3 

IF(AX. GT.  15.0)  NTEST  a  90.0  +  AX/2.0 

CBE00140  f 

CBE00150  t 

IF(N.LT. NTEST)  GO  TO  20 

CBE00160  N 

Is  /' 

WRITE(6,1001) 

CBE00170  ;•! 

K  :>: 

1001 

FORMAT(  IX,  'RANGE  OF  X  AND  N  IS  INCORRECT’) 

CBE00180 

r  '~ 

20 

RETURN 

CBE00190 

I  i 

N1 sN+1 

CBE00200  j* 

F  >* 

MA  a  AX  +  6.0 

CBE0021 0  W 

IF( AX. GE. 5 . 0 )MA  a  1.40»AX  +  60.0/AX 

CBE00220  ‘j 

J 

IX  a  AX 

CBE00230  fl 

\  ’» . 

MB  a  N+  IX/ 4  +  2 

CBE00240  J 

i  -:*‘ 

MO  a  MAXO(MA,MB) 

CBE00250  W 

• 

MMAX  a  NTEST 

CBE00260 

CBE00270  3 

DO  90  M  a  MO, MMAX, 3 

/  N- 

FM1  a  1. OB-28 

CBE00280  S 

✓ 

FM  a  0.0 

CBE00290  ■> 

J  \. 

ALPHA  a  O.DO 

CBE00300  £ 

i 

J1  a  1 

CBE0031 0  £ 

•t 

IF(M.  EQ . M/2*2)  J1  =  -1 

CBE00320  | 

P  ^ 

»  •> 

M2  a  M  -  2 

CBE00330  3 

DO  160  Ka1,M2 

CBE00340  Kl 

f  '• 

MK  a  M-K 

CBE00350 

/ 

XMK  a  MK  +  2 

CBEO0360  J 

5  >' 

_  •. 

BMK  a  XMK*FM1/X  -  FM 

CBE00370  cj 

FM  a  FM1 

cbeoo38o  r 

i 

FM1  a  BMK 

CBE00390  3 

>  ►> 

IF(MK-N-I.E).O.O)  CBJ  a  BMK 

CBE00400 

f  ^ 

J1  a  -J1 

CBE00410  d 

S  a  1  ♦  J1 

CBE00420  3 

I  <? 

j 

67 

i 

V' 


ALPHA  a  ALPHA  +  BMK*S 
BMK  a  2.0*FM1/X  -  FM 
IF( N.EQ.O)  CBJ  a  BMK 
ALPHA  a  ALPHA  +  BMK 
CBJ  a  CBJ/ ALPHA 

IF(CDABS(CBJ-CBPREV).LT.(CDABS(D*CBJ)))GO  TO  96 

CBPREV  a  CBJ 

RETURN 

END 

SUBROUTINE  CBJ(X,N,J) 

INPOTS  X,  N 

OOTPOT  J 

R  CO  TINE  TO  COMPOTE  COMPLEX  »16BESSEL  EUNCTION  JN(X) 

WHERE  X  IS  COMPLEX  *16 

THE  RCXJTINE  CBESSEL  (IB ESS. FTP  FROM  SOOTHAMPTON  IS  OSED 
IF  CABS(X)  IS  GREATER  THAN  7.0,  OTHERWISE  CBESJ  IS  OSED. 

ZERO  AND  FIRST  ORDER  BESSEL  FUNCTIONS  OF  COMPLEX  •16ARGOMENTS 
ARE  COMPUTED. 

N  a  0;  ZERO  ORDER  BESSEL  FUNCTION. 

N  a  1;  FIRST  ORDER  BESSEL  FUNCTION. 


SUBROUTINE  CBJ(X,N,J) 

IMPLICIT  REAL  »8  (A-H.O-Z) 

COMPLEX  »16  X, J, EMI ,EM, BMK, ALPHA, CBPREV 
COMPLEX  *16  Z,CHI,EZ,TO,T1,T2,T3,T4,T5,T6,T7,PO,QO 
REAL  *8  JR, JI 
AX  a  CDABS(X) 

IF(AX. GT.7.0)  GO  TO  100 
CBPREV  a  (0.0, 0.0) 

D  a  1.0E-04 

NTEST  a  20.0  +  10.0*AX  -  AX”2/3 
IF(AX. GT.15.0)  NTEST  =  90.0  +  AX/2.0 
IFCN.LT. NTEST)  GO  TO  20 
WRITE(6,1001) 

FORMAT( IX, 'RANGE  OF  X  AND  N  IS  INCORRECT’) 

RETURN 

N1  a  N  +  1 

MA  a  AX  ♦  6.0 

IF(AX. GE. 5.0)  MA  a  1.4D0*AX  +  60. DO/AX 
IX  a  AX 

MB  a  N  +  IX/4  +  2 

MO  a  MAXO(MA.MB) 

MMAX  a  NTEST 

DO  90  M  a  MO, MMAX, 3 

FM1  r  1.0E-28 

FM  a  0.0 

ALPHA  a  0.0 

JI  a  1 

IF(M.EH.M/2«2)  JI  a  -1 
M2  a  M  -  2 


CBE00430 
CBE00440 
CBE00450 
CBE0  0  460 
CBE00470 
CBE00480 
CBE00490 
CBE00500 
CBE0051 0 

TEM00010 

TEM00020 

TEM00030 

TEM00040 

TEM00050 

TEM00060 

TEM00070 

TEM00080 

TEM00090 

TEM00100 

TEM001 10 

TEMO 01 20 

TEM00130 

TEM00140 

TEM00150 

TEM00160 

TEM00170 

TEM00180 

TEM00190 

TEMO 0200 

TEM0021 0 

TEM00220 

TEMO 0230 

TEM00240 

TEM00250 

TEMO 0260 

TEM00270 

TEM00280 

TEMO 0290 

TEM00300 

TEMO 03 10 

TEM00320 

TEM00330 

TEM00340 

TEM00350 

TEMO 0360 

TEM00370 

TEM00380 

TEMO 03 90 

TEM00400 

TEM0041 0 

TEM00420 

TEM0  0  430 

TEM00440 

TEM00450 


y; 

to  . 

MVA'.'T' 

'"-  >.■: T^TTT  '.V.«.  ‘.'r  '.’Vv.'  '  ‘L’A'.VA*  <*V  'TVlr~-~<.Tvr».T'«.lr- 

iu 

?! 

t 

DO  160  K=1 , M2 

TEM00460 

1  /. 

HE  =  M  -  K 

TEM00470 

'*,» 

»'/. 

vV 

XMK  =  MK*2 

TEM00480 

i 

BMK  a  XMK*FM1/I  -  FM 

TEM00490 

4 

1 

FM  =  FM1 

TEM00500 

FM1  =  BMK 

TEM00510 

? 

IF(MK-N-1.H3.0)  J  =  BMK 

TEM00520 

\ 

*\ 

J1  a  -J1 

TEM00530 

:•: 

C-; 

S  a  1  +  J1 

TEM00540 

to 

160 

ALPHA  a  ALPHA  +  BMK*S 

TEM00550 

* 

BMK  a  2.0*FM1/X  -  FM 

TEM00560 

\ 

_-j 

0 

IF(N.EQ.O)  J  a  BMK 

TEM00570 

s 

r’ 

,N 

*• 

ALPHA  a  ALPHA  +  BMK 

TEM00580 

J  a  J/ALPHA 

TEM00590 

S 

to 

y. 

IF(CDABS( J-CBPHEV) .LT. (CDABS(D*J) ) )  GO  TO  200 

TEM00600 

XI 

90 

CBPREV  a  J 

TEM00610 

•3. 

V. 

GO  TO  200 

TEM00620 

V 

C 

BEGIN  CALCULATION  USING  LARGE  ARGUMENT  SERIES 

TEM00630 

C 

A.ANDS.EQU.9.2.5  (8  TERMS  ARE  USED) 

TEM00640 

•> 

100 

Z  a  X 

TEM00650 

EZ  a  8.0*Z> 

TEW)  0660 

IF(N.EQ.O)  GO  TO  189 

TEM00670 

f 

C 

TEM00680 

C 

FIRST  ORDER  CALCULATIONS 

TEM00690 

£>' 

.V 

C 

TEM00700 

CHI  a  Z  -  3.0*3.141 59265/4.0 

TEM00710 

% 

C 

TEM00720 

1 

C 

TERMS  IN  PO  AND  QO  SERIES  ARE  FORMED  BY 

TIM00730 

l 

C 

CHAIN  MULTIPLICATION 

TEM00740 

k. 

C 

TEM00750 

> 

TO  a  1.0 

TEM00760 

> 

•  • 

S. 

T1  a  T0*(3.0/EZ) 

TEW) 0770 

> 

V 

%*. 

T2  a  T1*(5.0/(2.0*EZ)) 

TEM00780 

m 

T3  a  T2*(21.0/(3.0*EZ)) 

TEM00790 

* 

1  *. 

T4  a  T3*(45.0/(4.0*EZ) ) 

TEM00800 

§ 

T5  a  T4*(77.0/(5.0*EZ) ) 

TEM00810 

■ 

T6  a  T5*( 117. 0/(6. 0*EZ)) 

TEM00820 

K 

T7  a  T6*( 165. 0/(7* D0*EZ)) 

TEM00830 

C 

TEM00840 

Cn 

C 

FORM  PO  a  1.0+T2-T4+T6 

TEM00850 

C 

TEM00860 

. 

PO  a  T6  -  T4 

TEM00870 

’.*• 

PO  a  PO  ♦  12 

TEM00880 

PO  a  PO  +T0 

TEM00890 

'.*■ 

C 

TEM00900 

»*• 

C 

FORM  QO  a  T1-T3+T5-T7 

TEM0091 0 

y 

c 

TEM00920 

f; 

QO  a  T5  -  T7 

TEM00930 

fc 

QO  a  QO  -  T3 

TEM00940 

* 

QO  a  QO  +  T1 

TEM00950 

/> 

GO  TO  199 

TEM00960 

£ 

c 

TEM00970 

j 

c 

ZERO  ORDER  CALCULATION 

TEM00980 

s 

c 

TEM00990 

s 

% 

189 

CHI  a  Z  -  3.14159265359/4.0 

TEM01000 

F 
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■i'.vv:*. 

r»  it-  V  '  v  *  v 

t 

C 

temoioio 

| 

c 

TERMS  IN  P0  AND  QO  SERIES  ARE  FORMED  BY 

TEM01020 

’.V 

c 

CHAIN  MULTIPLICATION. 

TEM01030 

ft 

c 

TEM0 1040 

\v 

•  V 

* 

TO  s  1.0 

TEM01050 

T1  a  TO/EZ 

TEM01060 

ft 

T2  a  T1 *(9.0/(2.0*EZ) ) 

TEM01070 

ft 

(fm 

T3  =  T2*(25.0/(3.0*EZ) ) 

TEM01080 

fa 

f: 

T4  a  T3f(49.0/(4.0»EZ)) 

TEM01090 

ft 

T5  a  T4«(81.0/(5.0«EZ)) 

TEM01100 

ft 

m 

L% 

T6  a  T5*( 121 .0/(6.0*EZ) ) 

TEMO 1110 

T7  a  T6  *C 169.0/(7.0*EZ) ) 

TEM01120 

iV 

c 

TEMO 1130 

< 

c 

FORM  PO  =  1.0  -  T2  +  T4  -  T6 

TEMO 11 40 

i 

ft 

c 

TEMO 1150 

ft. 

PO  =  T4  -  T6 

TEMO 1 160 

ft 

PO  =  PO  -  T2 

TEMO 1170 

ft 

PO  a  PO  +  TO 

TEMO 11 80 

ft 

IN 

c 

. 

TEMO 11 90 

c 

FORM  QO  a  -T1+T3-T5+T7 

TEMO 1200 

ft 

£• 

P 

c 

QO  a  T7  -  T5 

TEMO 1210 
TEMO 1220 

ft 

QO  a  QO  +  T3 

TEMO 1230 

E 

QO  a  QO  -  T1 

TEMO 1240 

ft 

199 

CONTINUE 

TEMO 1250 

ft 

ft 

J  a  CDSQRT( 2.0/ (3 . 1 41 59265359*Z) ) •( PO*CDCOS(CHI)- 

TEMO 1260 

2  Q0*CDSIN( CHI} ) 

TEMO 1270 

% 

C 

END  OF  CALCULATION  USING  SERIES  EXPRESSION. 

TEMO 1280 

K 

200 

CONTINUE 

TEMO 1290 

St 

XR  a  DREAL  (X) 

TEMO 1300 

y; 

ft 

XI  a  DIMAG(X) 

TEM01310 

ft 

JR  a  DREAL  (J) 

TEMO 1320 

ft 

JI  a  DIMAG (J) 

TEMO 1330 

& 

C 

FOR  ZEROTH  ORDER  FUNCTION : 

TEMO 1340 

9 

C 

-  IF  XRrO.O  OR  XI=0.0  THEN  J  IS  REAL  . 

TEMO 1350 

\V 

C 

-  THUS  SET  IMAG.  RESULT  TO  ZERO. 

TEMO 1360 

ft 

C 

IF( ( (XR.EQ.O.O) .OR. (XI.EQ.O.O) ) .  AND.  (N.EQ  .0) )  JI  a  0.0 

TEMO 1370 

TEMO 13 80 

& 

s.  , 

»_*  f 

C 

TEMO 1390 

& 

c 

FOR  FIRST  ORDER  FUNCTION 

TEMO 1400 

i  i 

77 

c 

-  IF  XI  a  0.0  THEN  J  IS  REAL 

TEMO 1410 

5c 

»  * 

*  • 

c 

-  THUS  SET  IMAG.  RESULT  TO  ZERO 

TEMO 1420 

w 

’  N1 

K.  * 

c 

TEMO 1430 

ft. 

c 

-  IF  XR  a  0.0  THEN  J  IS  IMAGINARY 

TEMO 1440 

*  V 

W 

V 

c 

-  THUS  SET  REAL  »8  RESULT  TO  ZERO 

TEMO 1450 

■**-  • 

t  ’ 

c 

TEMO 1460 

F 

IF((XI.EQ.0.0).AND.(N.EQ.1))  JI  a  0.0 

TEMO 1470 

pq 

,< 

.V 

A 

IF( (XR.EQ .0.0) .AND. (N. EQ .1 ) )  JR  a  0.0 

TEMO 1480 

J  a  DCMPLX( JR, JI) 

TEMO 1490 

RETURN 

TEMO 1500 

.yj 

VV 

END 

TEMO 1510 

d 

f 

• 

53 

£ 

p 

‘v 

ft 

V-' 

(-■ 

fe 

V-* 

►  , 

70 

• 
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APPENDIX  A. 2 


v:  c 

c 

\*  c 
c 
c 


c 


AMPL. FORTRAN  solves  for  the  amplitudes  of  the  wave  motion  developed  in 
the  section  Boundary  Value  Problem.  This  routine  calls  the  subroutines 
PROP,  SURIMP ,  LTRANR,  TRANSF  and  RTRANS  which  are  discussed  in  additional 
appendices.  AMPL  reads  the  following  input  from  a  data  file: 


VP1 

VP2 

VS1 

VS  2 

OMEGA 

THETO 

D 

RHOS 

PHI 

QSQ 

STATSF 

NTHER 

NVISC 

IARGT 

DM 

DE? 

D1 


p-wave  Velocity  in  Sand 
p-wave  Velocity  in  Clay 
s-wave  Velocity  in  Sand 
s-wave  Velocity  in  Clay 
Porosity 

Incidence  Angle  from  Normal 
Propagation  Depth 
Clay  Density 
Flow  Resistivity 
Tortuosity 

Dynamic  Shape  Factor  (s) 
Static  Shape  Factor  (n) 

Same  as  NTHER 

Set  to  $  for  this  work 

Not  used  for  IARGT  »  $ 

Layer  Thickness 
Propagation  Height 


For  each  frequency,  AMPL  calls  PROP  which  returns  the  fast  and  slow 
wave  propagation  constants,  then  the  constants  in  the  nine  boundary  condi¬ 
tion  equations  are  computed.  The  IMSL  routine  LEQTIC  is  called  and  returns 
the  wave  amplitudes.  Originally  AMPL  calculated  attenuation  coefficients, 
equation  73,  but  these  lines  (4710-4920)  are  no  longer  used. 


NUMERICAL  SOLUTION  FOR  WAVE  AMPL.  OF  A  BIOT  LAYER 


IMPLICIT  REAL  •  8  (A-H,0-Z) 

DIMENSION  WK(9) 

DIMENSION  ZREAL( 12) ,ZIMAG( 12) ,FRQ ( 12) , PHSGE0( 12) , PHSPRB( 12) 
REAL  *8  ATDBG{ 12) , ATDBP( 12) ,MAGG2( 12) ,MAGG1 (12) ,MAGP1 (12) 
REAL  »8  MAGP2( 12) ,MAGGE0( 12) ,MAGPRB( 12) 

COMPLEX  #16  HH1 ,CC1 ,MM1 
COMPLEX  *16  WA(99) 

DIMENSION  WA(9) 

COMPLEX  *16  TEMPI , TEMP2 , TEMP3 • TEMP4 , TEMP5 
COMPLEX  *16  Q3,Q4,Q5,C3,C4,C5 

COMPLEX  *16  J,L1 ,L2,L1S,L2S,C1 ,C2,S1 ,S2,Q1 ,Q2 
COMPLEX  *16  EXP1,EXP2,EXP3,EXP4,EXP5,EXP11,EXP12 


AMP00010 

AMP00020 

AMP00030 

AMP00040 

AMP00050 

AMP00060 

AMP00070 

AMP00080 

AMP00090 

AMP00100 

AMP00110 

AMP00120 

AMP00130 

AMP00140 

AMP00150 

AMP00160 

AMP00170 


ft 


COMPLEX  *16  EXP13,EXP14,EXP15 
COMPLEX  *16  VI , V2,Z1 ,Z2,M1 ,M2,F 
REAL  *8  OMEGA, KR,KF, KB BAR,K 
REAL  *8  MAGGO, MAGPB 
REAL  »8  MO , STATSF, NTHER, NV ISC, NPRIME 
REAL  *8  GB, GB2,LAM2 
COMPLEX  »16  A(9,9),B(9) 

REAL  *8  L0,L3,L4,L5,L0S,L3S,L4S,L5S,L0S0 

COMMON  /BLK3/  RHOF, GB, KBBAR , KR, RHOS, OMEGA, KF 
COMMON  /BLK2/  J,F,ETA,K,QS3 
COMMON  /BLK5/  70, CO, SO, Cl ,L0,L1 ,C2,L2,S3 ,L3 ,C3 
COMMON  /BLK4/  PHI 

COMMON  /BLK/  DM, IARGT, ICOMP, MO, STATSF, NTHER, NVISC, NPRIME 
COMMON  /BLK7/  HH1,CC1,MM1 
COMMON  /BLK8/  SI ,S2 
COMMON  /BLK8/  C4,S4,C5,L4,L5 
J=(O.DO, 1 .DO) 

EsDEXP(I.DO) 

TW0PI=4 . DO*DARSIN( 1 . DO ) 


DEP  IS  LAYER  THICKNESS 
D  IS  PROPAGATION  DEPTH 
D1  IS  A  PROPAGATION  HEIGHT 


WRITE(6,20) 

FORMATOX,  ’IMPDT  1  FOR  SARDIS,  2  FOR  ADDI  ACRES.') 

READ( 5 , * ) NTYPE 

IF( NTYPE  .GT.  2) GO  TO  9997 

READ( 2, 1 0 , END=99 ) VP1 , VP2 , VS1 , VS2 , OMEGA, THETO , D, RHOS, PHI,QS3 , 
2  , NTHER, NVISC, IARGT, DM, DEP, D1 
IF( NTYPE  .NE.  1)G0  TO  409 

READ( 1 , 1 0, END=409)VP1 , VP2 , VS1 , VS2 , OMEGA, THETO ,D, RHOS, PHI,  QS3 
2  STATSF, NTHER, NVISC, IARGT, DM, DEP, D1 

FORMAT(/,1X,7(E10.3,1X),///,1X,6(E10.3f1X),I3,//,1X,3(E10.3 
WRITE ( 1 1 ,21 ) VP1 , VP2 ,VS1 , VS2 , OMEGA, THETO , D, PHI, Q S3 , STATSF, 

2  NTHER, NVISC, IARGT,  EM, D1, DEP 

FORMATOX,  'VP1  s  »,E10.3,2X, 'VP2  =  » ,E1 0.3 ,2X, ' VS1  =  ' 

2  ,E10.3,3X, ' VS2  s  »,E10.3,/,1X, 'OMEGA  s  ',E10.3,3X, 

2  'THETO  =  *,E10.3,3X, '+Z  =  ' ,E10.3 ,/,1X, 'ELCW  RESIST  =• 

2  , El 0.3 ,3X, 'TORTUOSITY  s  ' ,E10. 3, 3X, 'STATSF  =  ' 

2  ,E10.3,/,1X, 'NTHER  s  ', El 0.3, 3X, 'NVISC  =  ', 

2  E10.3.3X, 'IARGT  =  »,I2,3X,'DM=  ' , El  0.3 ,3X,/, IX, »-Z  =  »,E1 

2  , 'LAYER  THICKNESS  =  '.E10.3) 

FR  IS  FREQUENCY 
THETO  IS  INCIDENT  ANGLE 

DO  9  Is  1,300 
FRQ  ( I ) = I 
CONTINUE 
FRslO.DO 


AMP00180 
AMP00190 
AMP00200 
AMP0021 0 
AMP00220 
AMP00230 
AMP00240 
AMP00250 
AMP00260 
AMP00270 
AMP00280 
AMP00290 
AMP00300 
AMP003 1 0 
AMP00320 
AMP00330 
AMP00340 
AMP00350 
AMP00360 
AMP00370 
AMP00380 
AMP00390 
AMP00400 
AMP0041 0 
AMP00420 
••••••AMP00430 

AMP00440 
AMP00450 
AMP00460 
AMP00470 
AMP00480 
STATSFAMP00490 
AMP00500 
AMP0051 0 
,  AMP00520 

AMP00530 
,1X) )  AMP00540 
AMP00550 
AMP00560 
AMP00570 
AMP00580 
AMP00590 
AMP006 00 
AMP00610 
0.3,3XAMP00620 
AMP00630 
••••••AMP00640 

AMP00650 

AMP00660 

AMP00670 

AMP00680 

AMP00690 

AMP00700 

AMP00710 


o  o  o  n  o  o  o  n  n  o  n  o  o 


C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


LTHETI=THETO 

THETO = TWOPI«THETO / 36 0 . DO 
D3=D 


D  IS  PROPAGATION  DEPTH  BELOW  SURFACE. 
D1  IS  PROPAGATION  DEPTH  ABOVE  SURFACE/ 


WAVE  NUMBERS 


0  INCIDENT  WAVE  IN  AIR  AT  THETO. 

1  FAST  WAVE  IN  POROUS  MEDIUM 

2  SLOW  WAVE  IN  POROUS  MEDIUM 

3  SHEAR  WAVE  IN  POROUS  MEDIUM 

4  SHEAR  WAVE  IN  3RD  MEDIUM 

5  FAST  WAVE  IN  3RD  MEDIUM 


imHIHHHIIIHHHIlHHtHHHHCHHHIHi 

VO  TO  V5  ARE  THE  WAVE  VELOCITIES  IN  CGS  UNITS 

VO=  SOUND  SPEED  IN  AIR 

VI =  BIOT  FAST  WAVE  VELOCITY 

V2=  BIOT  SLOW  WAVE  VELOCITY 

V3=  SHEAR  WAVE  VELOCITY  IN  LAYER 

V4=  SHEAR  WAVE  VELOCITY  IN  LOWER  MIDIUM 

V5=  COMPR ESSIONAL  WAVE  VELOCITY  IN  LOWER  MEDIUM 


DENSITIES  OF  MATERIALS  FOR  THE  SYSTEM: 
DENSITY  OF  AIR,  RHOF  (CGS  UNITS) 
RHOF=1.204D-03 


DENSITY  OF  CLAY,  RHOS  AND  MEAS  VAL.  =  2.65  GCM*«-3 


DENSITY  OF  POROUS  LAYER,  RHO 


RHO=RHOS+OMEGA*( RHOF -RHOS) 


GB  AND  GB2  ARE  THE  SHEAR  MODULI  FOR  THE  LATER  AND  CLAY 

GB  =  VS1  #VS1  *RHO 
C 

GB2  »  VS2*VS2*RH0S 


AMP00720 
AMP00730 
AMP00740 
AMP00750 
AMP00760 
AMP00770 
AMP00780 
AMP00790 
AMP00800 
AMP0081 0 
AMP00820 
AMP00830 
AMP00840 
AMP00850 
AMP00860 
AMP00870 
AMP00880 
AMP00890 
AMP00900 
AMP00910 
AMP00920 
AMP00930 
AMP00940 
AMP00950 
AMP00960 
AMP00970 
AMP00980 
AMP00990 
AMP01000 
AMP01010 
AMP01020 
AMP01030 
AMP01040 
AMP01050 
AMP01060 
AMP01070 
AMP01080 
AMP01090 
AMP01100 
AMP01110 
AMP01120 
AMP01130 
AMP01140 
AMP01150 
AMP01160 
AMP01170 
AMP01180 
AMP01190 
AMP01200 
AMP01210 
AMP01220 
AMP01230 
AMP01240 


73 


V  ' 


f5 


c. 


c 

c 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cc 

30 

c 

c* 

c 

c 

11 

c 

c 

c 

c 


c 

50 


c 

c 

c 

c 

c 

c 

65 


KBBAR  =  VP1»VP1»RH0  -  4.D0/3.D0«GB 
KR  =  VP2*VP2*RH0S  -  4.D0/3.D0*GB2 


EQUATIONS  TO  CALCULATE  H,C,M  FOR  THE  POROUS  MEDIUM 
ARE  TAKEN  FROM  STEVE  ADDISONS  SEMINAR. 

H  IS  HH1 ,C  IS  CC1,  M  IS  MM1 


FOR  THE  LATER:  BULK  MOD  FOR  THE  GRAINS  IS  SAME  AS  CLAY,  KR 
BULK  MOD  FOR  THE  FRAME  IS  INPUT,  KBBAR 
BULK  MOD  OF  FLUID  =ATM.  PRES.  1.01E6  KF 

KF=1 .01D06 


TEMP=KR*( 1 .D0+0MEGA»(KR/KF-1 .DO) ) 

HH1 = (  ER- KBBAR ) «»2 . DO/ ( TEMP- KBBAR ) +KBB AR+GB*4 .DO/3. DO 
CC1=KR*(KR-KBBAR)/( TEMP- KBBAR ) 

MM1 =KR *»2 . DO/ ( TEMP-KBBAR ) 


NRITE(6,30)HH1 ,CC1,MM1 

FORMATdX,  'H=*  ,2E13.6,2X,  'C='  ,2E13.6,2X,  'M='  ,2E13-6) 


DO  101  J1 =5,10,5 
DO  201  1=1,30 
CALL  PR0P(FR,L2,L1) 


WRITE(6,30)HH1 ,CC1,MM1 

VI =TW0PI«FR/L1 

V2=TW0PI«FR/L2 

WRITE( 11 ,50)  L2,L1 

F0RMAT(5X, 'BIOT  WAVE  CONSTANTS’,/ 

2  ,1  OX, 'SLOW  WAVE  CON.  =  ’  ,E1 3.6 , ' , '  ,E1 3 .6 , 

2  /,10X, 'FAST  WAVE  CON,  =  ' ,E1 3 .6 , ' , ' ,E1 3 .6 ,/) 

V0=3-43D4 


WRITE( 6 , 65 ) KBBAR , KR , GB, GB2 

F0RMAT( IX, ’KBBARs’ ,E13.2,2X, ’KR=’ ,E13.2,2X, 'GB=' , El  3. 2, 2 X 
2  , *GB2s' , El  3.2,///) 


AMP01250 
AMP01260 
AMP0 1 27  0 
AMP01280 
AMP01290 
AMP01300 
AMP01310 
AMP01320 
AMP01330 
AMP01340 
AMP01350 
AMP01360 
AMP01370 
AMP01380 
AMP01390 
AMP01400 
AMP01410 
AMP01420 
AMP01430 
AMP01440 
AMP01450 
AMP01460 
AMP01470 
AMP01480 
AMP01490 
AMP01500 
AMP01510 
AMP01520 
AMP01530 
AMP01540 
AMP01550 
AMP01560 
AMP01570 
AMP01580 
AMP01590 
AMP01600 
AMP01610 
AMP01620 
AMP01630 
AMP01640 
AMP01650 
AMP01660 
AMP01670 
AMP01680 
AMP01690 
AMP01700 
AMP01710 
AMP01720 
AMP01730 
AMP01740 
AMP01750 
AMP01760 
AMP01770 
AMP01780 


V 


L 


c 

AMP01790 

if; 

c 

AMP01800 

■V» 

c»««* 

'••••AMP01810 

c 

AMP01 820 

fl 

c 

Ml  AND  M2  ARE  THE  RATIOS  OF  THE  RELATIVE 

WAVE  TO 

THE 

AMP01830 

f? 

c 

MATRIX  WAVE(FAST,  SLOW) 

AMP01840 

c 

AMP01850 

c 

EQUATIONS  FOR  Ml  AND  M2  TAKEN  FROM  GEERTSMA  AND  SMIT 

AMP0 1 86  0 

A’. 

c 

GEOPHYSICS  V.26  APR  1961 

AMP01870 

c 

ASPECTS  OG  ELASTIC  WAVE  PROG.  IN  ELD.  SAT. 

POROUS 

SOL  IDS  • 

AMP01880 

c 

AMP01890 

*4 

V'„ 

C**H* 

issssamP01900 

“\ 

c 

AMP01910 

21=  HH1 *L1 •L1/(RH0*TW0PI#TW0PIiFR*FR) 

AMP01920 

\/ 

c 

AMP01930 

vV 

Z2=  HH1 *L2*L2/(RH0BTW0PI#TW0PIiFRtFR) 

AMP01940 

c 

AMP01950 

>  •' 

M1=  (HH1 *(Z1-1 .DO) )/(Z1 *CC1-RH0F*HH1/RH0) 

AMP01960 

iV 

c 

AMP01970 

M2=  (HH1 »(Z2-1 .DO) )/(Z2#CC1-RH0F*HH1 /RHO) 

AMP01980 

c 

AMP01990 

•.- 

c 

WRITE( 11 , 40)M1 ,M2 

AMP02000 

40 

F0RMAT(1X,  »M1  =  \E13.6,E13.6,2X,  'M2=* ,  El  3-6 

,E13.6,///) 

AMP02010 

c 

AMP02020 

,V 

>*** 

AMP02030 

V, 

c 

AMP02040 

C  LAME 

'  S  CONSTANTS  FOR  THE  LOWER  MEDIUM  ARE  LAM2, 

GB2 

AMP02050 

Iftf 

C 

AMP02060 

E 

LAM2=KR  +  ( 4 . DO/3 • DO )  *GB2 

AMP02070 

C 

AMP02080 

C 

AMP02090 

V* 

C 

AMP021 00 

V 

C 

AMP02110 

c*«« 

AMP02120 

• 

1  c 

AMP02130 

(*  " 

c 

LO  -  L5  ARE  THE  WAVE  NUMBERS 

AMP021 40 

c 

AMP021 50 

.*, 

c 

AMP02160 

LO=TWOPI*FR/VO 

AMP02170 

L3=TW0PI*FR/VS1 

AMP021 80 

1  *• 

L4=TW0PI«FR/VS2 

AMP02190 

»  ', 

L5 = TWOPI *FR/ VP2 

AMP02200 

h  * 

c 

AMP0221 0 

c 

AMP02220 

/• 

LOS=LO*LO 

AMP0223  0 

i  ' 

LI S=L1 *L1 

AMP02240 

L2S=L2*L2 

AMP02250 

. " 

> 

L3S=L3#L3 

AMP02260 

■i 

f 

L4S=L4*L4 

AMP02270 

L5S=L5*L5 

AMP02280 

o  C 

IS 

AMP02290 

t  * 

c 

AMP02300 

* 

c 

SO, SI, CO, Cl,  ETC  ATE  THE  SINES  AND  COSINES 

OF  THE 

AMP02310' 

„ 

c 

INCIDENT,  REELECTED  AND  REFRACTED  ANGLES. 

AMP02320 

*.*.  c 

*  V 


ii 


75 


I; 


c 


1 


'r  ^ 


£ 


V 

v 


c 

c 

c 

c 

c 

c 

c 


i 


« 


c 

c 

c 

c 

c 

c 

c 


AMP02330 

AMP02340 

SOsDSIN(THETO) 

AMP02350 

COsDCOS(THETO) 

AMP02360 

S3Q=S0*S0 

AMP02370 

Cl =CDS3  RT( 1 . DO-LOS*SSQ/L1  S) 

AMP023  80 

C2=CDSQRT(  1  .DO-LOS*SSQ/L2S) 

AMP02390 

TEMP3=1 . DO-LOS*SSQ/L3S 

AMP02400 

TEWP4=1 .D0-L0S*SSQ/L4S 

AMP0241 0 

TEMP5=1 .D0-L0S»SSQ/L5S 

AMP02420 

C3=CDSQRT(TEMP3) 

AMP02430 

C4  sCDSQ  RT ( TEMP4 ) 

AMP02440 

C5=CDSQ RT(TEMP5 ) 

AMP02450 

LOSO=LO*SO 

AMP02460 

SI =LO SO/LI 

AMP02470 

S2=L0S0/L2 

AMP02480 

S3=LOSO/L3 

AMP02490 

S4=L0S0/L4 

AMP02500 

S5=L0S0/L5 

AMP02510 

AMP02520 

AMP02530 

Q1  -  Q5  IS  THE  HAVE  NUMBER  •  COS( ANGLE)  •  LATER 

THICKNESS 

AMP02540 

AMP02550 

DEP  IS  THE  DEPTH  OF  LAYER 

AMP02560 

AMP02570 

Q1=L1*C1*DEP 

AMP02580 

AMP02590 

Q2=L2*C2*DEP 

AMP026  00 

Q3sL3#C3«DEP 

AMP02610 

Q4=L4«C4«DEP 

AMP02620 

Q5=L5*C5*DEP 

AMP02630 

WRITE(6,*)Q1 ,Q2,Q3jQ4,Q5 

AMP02640 

AMP026  50 
AMP02660 

EXPI  OR  EXPII  ARE  THE  EXPONENTIALS  IN  THE  B.  C. 

EQUATIONS 

AHP02670 

AMP02680 

EXP 11  =CDEXP(-J*Q1 j 

AMP02690 
AMP027  00 
AMP02710 

EXPI =CDEXP( J*Q  1 ) 

AMP02720 

£ 


I  4 

t’- 


<\ 

r.' 


c 
c 

V'  C 

t  % 


EXP 12=CDEXP(-J*Q2) 
EXP2=CDEXP( J*Q2) 
EXP3=CDEXP(J«Q3) 
EXP 1 3 =CDEXP ( - J "Q  3 ) 
EXP4sCDEXP(J«Q4) 
EXP14=CDEXP(-J*Q4) 
EXF5=CDEXP( J«Q5) 
EXP15=CDEXP(-J«Q5) 


THE  A  SUB  II' S  ARE  THE  COEF.  IN  THE  B.  C.  EQUATIONS. 


AMP02730 
AMP02740 
AMP02750 
AMP02760 
AMP02770 
AMP02780 
AMP02790 
AMP02800 
AMP0281 0 
AMP02820 
AMP02830 
AMP02840 
AMP02850 
AMP02860 


\Vj| 

J 

t'4 

m 

tf>V 

VlSi 

S3 

W! 


I 


rwzr-'r~w.'v-*r 

> 

TWTTO*  V  V-  ■  i  "l  "ST*  i.WlH’W-1 

nu nuft 

r'  l  ■ 

*  t-* 

J  L 

» 

s 

% 

i  rt 

A( 1 , 1 )=( 1 .DO-MI )*C1 

AMP02870 

r.  F.n 

A(1 ,2)=-A(  1,1 ) 

AHP02880 

J 

N 

N 

A(  1 »3)  =  ( 1 .D0-M2)*C2 

AMPQ2890 

ri 

c 

ft  f 

A(1,4)=-A(1,3) 

AMP02900 

J 

I  R 

A(1,5)=S3 

AMP0291 0 

A(1,6)=-S3 

AMP02920 

9 

it  ■ 

A( 1 ,7) =( 0.D0,0 . DO) 

AMP02930 

¥ 

A( 1 ,8)=(0.D0,0.D0) 

AMP02940 

.*  *■  *  • 

A(  1 ,9)SC0 

AMP02950 

J 

S  f» 

A(2, 1 )=2.D0*L1 *S1 «C1 

AHP02960 

L 

A(2,2)=-A(2, 1 ) 

AMP02970 

'*  i/ 

A( 2 , 3 ) =2 . D0*L2  *C2*S2 

AHP02980 

V 

< 

< 

A(2,4)=-A(2,3) 

AMP02990 

3 

>  $ 

A( 2 , 5 ) =- ( L3»C3«C3-L3«S3  #S3 ) 

AMP03000 

> 

'  & 

A(2,6)=-A(2,5) 

AHP03010 

s 

| 

A(2,7)=(0.DO,0 .DO) 

AMP03020 

> 

A(2,8)=(0.D0,0.D0) 

AMP03030 

; 

j  >; 

A(2,9)=(0.D0,0.D0) 

AMP03040 

A(3,1)=L1 •(HH1-M1«CC1-2.D0«GB»S1  «S1 ) 

AMP03050 

z 

•4I 

-A.  ,  - , 

A(3,2)=A(3,1) 

AMP03060 

j 

£  £ 

A(3,3)=L2*(HH1-M2»CC1-2.D0«GB*S2»S2) 

AMP03070 

i  f 

A(3,4)=A(3,3) 

AMP03080 

F 

\ 

A(3,5)=2.D0*GB«L3*S3*C3 

AMP03090 

*  * 

-*  .  V 

A(3»6)=A(3»5) 

AMP03100 

>  /- 
.*  *  *-f 

A(3,7)=(O.DO,O.DO) 

AMP03110 

► 

l] 

• 

\ 

A(3»8)=(0.D0,0.D0) 

AMP03120 

> 

A(3,9)=-KF*L0 

AMP03130 

1  I 

A(4,1 )=L1*(CC1-M1*MM1) 

AMP03140 

L 

m  V 

A(4,2)=  A(4,1 ) 

AHP03150 

«* 

-  _  • 

A(4,3)=L2*(CC1-M2»MM1) 

AMP03160 

r. 

* ,  > 

■*  «* 

A(4,4)=A(4,3) 

AMP03170 

i 

•rf  -  ■• 

/ 

A(4,5)=(O.DO,O.DO) 

AMP03 1 80 

t 

i  _ 

A(4,6)=(0.D0,0.D0) 

AMP03190 

1  f 

A(4,7)=(0.D0,0.D0) 

AMP03200 

S 

c  t-. 

A(4,8)=(0.D0,0.D0) 

AMP03210 

y 

s 

% 

A(4,9)=-KF*L0 

AMP03220 

Hi 

A(5,1)=C1*EXP1 

AMP03230 

Q  '•:* 

A(5»2)=-C1*EXP1 1 

AMP03240 

5 

g 

A(5,3)=C2»EXP2 

AMP03250 

3 

A(5,4)=-C2«EXP12 

AMP03260 

A(5,5)=S3»EXP3 

AMP03270 

F 

[•  ■'■ 

A(5,6)=-S3»EXP13 

AMP03280 

fc 

r, 

A(5»7)=-C5*EXP5 

AMP03290 

■;  ‘X 

A(5,8)=-S4*EXP4 

AMP03300 

/j 

A 

i  <* 

A(5,9)=(O.DO,O.DO) 

AMP03310 

K 

r 

A(6,1 )=S1 »EXP1 

AMP03320 

a 

> 

A(6,2)=S1 •EXPII 

AMP03330 

&j 

i  -  •; 

A(6 ,3) =S2*EXP2 

AMP03340 

v» 

A(6»4)sS2*EXP12 

AMP03350 

s 

/ 

A(6,5)s-C3tEXP3 

AMP03360 

5 

«  y 

A(6,6)=-C3*EXP13 

AMP03370 

L 

W 

A(6,7)=-S5»EXP5 

AMP03380 

% 

^  (V 

A(6 ,8)sC4*EXP4 

AMP03390 

s 

A(6,9)s(O.DO,O.DO) 

AMP03400 

V 

5  >:■ 

A(7» 1 ) » A( 3 » 1 )*EXP1 

AHP03410 

t  ^ 

ri!>V*V 


n  t1 


o  o  o  o  o  o  o  o  -*  n  o  oonoono 


AMP0342G 

AMP03430 

AMP03440 

AMP03450 

AHP03460 

AMP03470 

AMP03480 

AMP03490 

AMP03500 

AMP03510 

AMP03520 

AMP03530 

AMP03540 

AMP03550 

AMP03560 

AMP03570 

AHP03580 

AMP03590 

AMP03600 

AMP03610 

AMP03620 

AMP03630 

AMP03640 

AMP03650 

AMP03660 

AMP03670 

AMP03680 

'AMP03690 

AMP03700 


WHITING  AMPLITUDE  RATIOS  AMP03710 

AMP03720 

DO  280  Is 1,9  AMP03730 

WRITE(6,1 10) (A(I, K1 ) ,K1=1 ,9)  AMP03740 

280  CONTINUE  AMP03750 

AMP03760 

10  FORMAT( 1X,9(E1 8.8, ,,,,E14.8,/))  AMP03770 

AMP03780 

**************************************************************  AMP03790 

AMP03800 

AMP03810 

B( 1 )=C0  AMP03820 

B(2)s(0.D0,0.D0)  AMP03830 

B(3)=KF*L0  AMP03  840 

B(4)=KF*L0  AMP03850 

B(5)=(0.D0,0.D0)  AMP03860 

B(6)=(0.D0,0.D0)  AMP03870 

B(7)=(0.D0,0.D0)  AMP03880 

B( 8)3(0. DO, 0. DO)  AMP03890 

B( 9)=(0.D0,0 .DO)  AMP03900 

WRITE(6,*)(B(I),Is1,9)  AMP03910 

AMP03920 

AMP03930 


AMP03940 


A(7,2)=A(3,2)*EXP11 

A(7,3)sA(3,3)#EXP2 

A(7»4)=A(3»4)*EXP12 

A(7 , 5 ) =2 . D0*GB*L3  *S3*C3*EXP3 

A<7 , 6 ) s2 . D0*GB*L3  *S3  * C3*EXP 1 3 

A( 7 17) =-L5  *EXP5* ( LAM2-2 . D0*GB2«C5  *C5 ) 

A( 7 , 8 ) s2 . D0«GB2«L4*S4 «C4 «EXP4 

A(7,9)=(O.DO,O.DO) 

A(8,1)s2.D0*GB*L1*S1*C1*EXP1 

A( 8 , 2 ) =-2 . D0*GB*L1 *S1 *C1*EXP11 

A( 8, 3 ) =2 . D0*GB*L2*S2*C2*EXP2 

A(8,4)s-2.D0*GB*L2*S2*C2*EXP12 

A(8,5)s-GB*L3*(C3*C3-S3*S3)*EXP3 

A( 8,6) =GB*L3  *( C3*C3-S3  *S3 ) *EXP 1 3 

A(8f7)=-2.D0*GB2*L5*C5*S5*EXP5 

A( 8, 8 ) sGB2*L4*( C4«C4-S4«S4 ) *EXP4 

A(  8,9)  =  (0.D0,0 .DO) 

A( 9 » 1 ) =M1 *C1 *EXP1 
A(9»2)=-M1  *C1*EXP11 
A(9,3)=M2*C2*EXP2 
A( 9 , 4) s-M2*C2*EXP 1 2 
A(9i5)=(O.DO,O.DO) 

A(9, 6)3(0.00,0. DO) 

A(9,7)=(0.D0,0.D0) 

A(9, 8)3(0.00,0. DO) 

A(9»9)=(0.D0,0 .DO) 


C  AMP03950 

C  SOB ROD TINE  LEQ  TIC  AHP03960 

C  AMP03970 

C  SOLVES  SIMULTANEOUS  ALGREBRAIC  EQUATIONC  AMP03980 

C  AMP03990 

C  AMP04000 

c  •iH»*i»H«*»im»***f»*H«i*ii»»»imfi*ii»m»»ii»*HHH««  AMP04010 

IA=9  AMP04020 

IB=9  AMP04030 

N=9  AMP04040 

UOB=0  AMP04050 

M=1  AMP04060 

CALL  LEQ2C(A,  N,  IA,  Bt M,  IB,  UOBt WA, WK,  IER)  AMP04070 

C  AMP04080 

C  AMP04090 

C  MIHHMMlfmHIHHIIHHUHHIHIHHHIIIHIIHHHMH  AMP041 00 

C  AMP04110 

C  AMP041 20 

IFRHJ=FR  AMP04130 

C204  WRITE(1 1,205)  AMP04140 

205  FORMAT(5X, 'FREQUENCY' ,5X' INCIDENT  ANGLE*)  AMP04150 

C  WRITE (1 1,21 0)IFREQ,LTHETI  AMP04160 

210  FORMAT(6X, 14, 15X, 12,//)  AMP04170 

C  WRITE( 1 1 ,300)  AMP04180 

300  FORMAT(9X,  'FAST  WAVE  (DOWN) • ,20X, *  FAST  WAVE  (UP)*)  AMP04190 

C  WRITE( 1 1 ,305)B( 1 ) ,B( 2)  AMP04200 

305  FORMAT( IX, El  4. 8, IX, ' , * ,1X, El  4.8, 5X, El  4.8, IX, ' , ' ,1X,E1 4.8,//)  AMP04210 

C  WRITE ( 1 1 ,310)  AMP04220 

310  FORMAT(9X, 'SLOW  WAVE  (DOWN)*, 21X, *SLCW  WAVE  (UP)’)  AMP04230 

C  WRITE( 1 1 ,305)B( 3) ,B( 4)  AMP04240 

C  WRITE( 11 ,320)  AMP04250 

320  FORMAT(9X, 'SHEAR  WAVE  (DOWN)', 20X, 'SHEAR  WAVE  (UP)')  AMP04260 

C  WRITE(11,305)B(5),B(6)  AMP04270 

C  WRITE(1 1,330)  AMP04280 

330  FORMAT(9X, 'TRANS.  COMP.  WAVE’ ,20X, 'TRANS.  SHEAR  WAVE' )  AMP04290 

C  WRITE(11,305)B(7),B(8)  AMP04300 

WRITE(6,340)  AMP04310 

340  F0RMAT( 9X, ' REFL  WAVE')  AMP04320 

WRITE(6,305)B(9)  AMP04330 

C  WRITE(6,»)(B(I),I=1,9)  AMP04340 

C  AMP04350 

Cim«m«H«HmmiHi««mHHHi«i«imHHHiiimimH«fmmHAHP04360 


CALL  SUFIMP( B,  LO , KF, RHOF, VO , CO , FR ,  ZREAL 1 ,ZIMAG1) 

AMP04370 

ZREAL(I)sZREAL1 

AMP04380 

ZIMAG(I) sZIMAGI 

AMP04390 

c 

IF(D.GT.DEP)  GO  TO  1090 

AMP04400 

c 

CALL  TRANSF( B,M1,M2,FR,D,D1, MAGGO , MAG PB ,  PHSGO ,  PHSPB ) 

AMP04410 

c 

CALL  RTRANS(B, FE,D, D1 , RMAG, RPHAS) 

AMP04420 

GO  TO  1095 

AMP04430 

Cl  090 

CALL  LTRANF( B, LO , RHOF, KF, FR, DPRIME , MAGGO, PHSGO) 

AMP04440 

1095 

CONTINUE 

AMP04450 

MAGGEO( I ) = MAGGO 

AMP04460 

C 

WRITE  ( 6 , 1 3 )  FR ,  MAGGO ,  MAG  PB 

AMP04470 

Cl  3 

F0RMAT(1X,3(E14.3)) 

AMP04480 

i 


MAGPRB* I) aMAGFB 

AMP04490 

PHSGE0(I)=PHSG0 

AMP04500 

PHSPRB(I)=PHSFB 

AMP04510 

c 

AMP04520 

c 

WRITING  VALUES  FOR  TRANSFER  FUNCTION  AND  SURFACE  IMPEDANCE. 

AMP04530 

c 

AMP04540 

c 

WRITE*  6 , 256 ) MAGGEO* I ) , PHSGEO* I ) , MAGPRB ( I ) , PHSPRB ( I ) , Z  REAL ( I ) 

AMP04550 

c 

2 

,ZIMAG(I) , RMAG, RPHAS 

AMP04560 

256 

FORMAT* IX, 8 El  2. 4) 

AMP04570 

WRITE* 6 , 257 ) ZREAL * I) , ZIMAG ( I ) 

AMP04580 

257 

FORMAT* IX, 2E1 2.4) 

AMP04590 

505 

FORMAT* IX, ‘MAD  FOR  VERT  GEO.  =  '.E14.8,/ 

AMP04600 

2 

,1X, ‘PHAS  FOR  VERT  GEO.  =  »,E14.8,// 

AMP046 10 

2 

,1X, 'MAG  FOR  PROBE  =  »,E14.8,/, 

AMP04620 

2 

IX, 'PH AS  FOR  PROBE  =  »,E14.8,// 

AMP04630 

2 

, IX, 'REAL  PART  OF  IMPEDANCE  =  *,E14.8,/, 

AMP04640 

2 

IX, *IMAG  PART  OF  IMPEDANCE  =  '.E14.8,//) 

AMP04650 

C 

AMP04660 

C 

MAGG2 ( I ) =MAGGE0 ( I ) 

AMP04670 

C 

MAGP2(I)=MAGPRB(I) 

AMP046  80 

FR  a  FR  +  10.0D0 

AMP04690 

201 

CONTINUE 

AMP04700 

C 

IF* J1 . GE. 10)GO  TO  1100 

AMP04710 

C 

DO  301  11=1,12 

AMP04720 

C 

MAGG1 ( 11 )=MAGG2( 11 ) 

AMP04730 

C 

MAGPKI1  )=MAGP2(I1 ) 

AMP04740 

C301 

CONTINUE 

AMP04750 

C 

D  =  D  +  1 .DO 

AMP04760 

C 

FR  a  25. DO 

AMP04770 

C101 

CONTINUE 

AMP04780 

1100 

CONTINUE 

AMP04790 

C 

FR  =  25. 

AMP04800 

C 

DELD  =  D  -  D3 

AMP0481 0 

C 

DO  401  1=1,12 

AMP04820 

C 

ATDBG(I)=*-1.D0)*(DL0G(MAGG2(I))-DL0G*MAGG1*I)))*20.D0»DL0G10(E)/AMP04830 

C 

2 

DELD 

AMP04840 

C 

ATDBP(I) =*-1 .DO) *(DL0G(MAGP2(I) )-DL0G*MAGP1  *1) ) ) *20.DO#DLOG1 0(E)/AMP04850 

C 

2 

DELD 

AMP04860 

c 

FR  a  FR  +  25. DO 

AMP04870 

C401 

CONTINUE 

AMP04880 

C 

WRITE* 11, 203) 

AMP04890 

203 

FORMAT* IX, ’FREQ,  ATTENUATION  IN  DB  OF  PROBE  AND  VERT.  GEO.') 

AMP04900 

C 

WRITE* 11 ,506) (FRQ(I) , ATDBP(I) , ATDBG(I) , 1=1 ,12) 

AMP04910 

506 

FORMAT* IX, 3(E1 4. 8, IX)) 

AMP04920 

9997 

STOP 

AMP04930 

END 

AMP04940 

% 

£ 

C 
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APPENDIX  A. 3 


The  routine  SUF IMP. FORTRAN  computes  the  real  (ZREAL)  and  imaginary 
(ZIMAG)  parts  of  the  acoustic  surface  impedance  relative  to  air  using 
equation  67.  The  input  arguments  passed  are  all  defined  in  AMPL. FORTRAN 
(see  Appendix  A. 2). 


SUBROUTINE  TO  CALCULATE  THE 
REAL  AND  IMAGINARY  PARTS  OF 
THE  SURFACE  IMPEDANCE. 

SUBRGU TINE  SJ FIMP(B, LO , KF,  RHOF, VO , CO , FR , ZREAL , ZIMAG ) 
IMPLICIT  REAL  *8  (A-H,0-Z) 

COMPLEX  *16  B(9) 

COMPLEX  *16  J , PRES, VEL , ZNOT 
REAL  *8  KF.LO 

CALCULATE  PRESSURE  AND  VELOCITY  AT  THE  SURFACE. 
J=(0.D0,1.D0) 

TWOPI=  4 • DO*DARSIN ( 1 .DO) 

PRES  a  RHOF*VO*CO*( 1 .DO+BC 9 ) ) *(-J*TWOPI*FR) 

PRES  a  - J *KF*LO  * ( 1 . D0+B( 9 ) ) 

VEL=C0*( 1 .D0-B( 9) )*(-J*TWOPI*FR) 

ZNOTs  PRES/VEL 

ZREALaDREAL ( ZNOT) / ( RHOF*VO ) 

ZIMAGaDIMAG ( ZNOT) /( RHOF*VO ) 

WRITE ( 6 , 1 3 ) ZREAL , Z IMAG 
FORMAT( IX, ' IMPEDANCE  a  «,2(E14.3)) 

RETURN 

pirn 


SUFOOOIO 

SUF00020 

SUF00030 

SUF00040 

SUF00050 

SUF00060 

SUF00070 

SUF00080 

SUF00090 

SUF00100 

SUFOO1 10 

SUF00120 

SUF00130 

SUF00140 

SUF00150 

SUF00160 

SUF00170 

SUF00180 

SUF00190 

SUF00200 

SUF00210 


r_\ 


APPENDIX  A. 4 


The  subroutine  TRAN SF. FORTRAN  computes  the  acoustic  and  vertical  seis¬ 
mic  transfer  functions  for  the  porous  sand  using  equations  77  and  79, 
respectively.  TRANSF  returns  the  magnitudes  (DMAG,  DMMAG)  and  phases 
(PHAS,  PHAS2)  of  these  transfer  functions.  Input  parameters  are  defined  in 
AMPL. FORTRAN. 


LET  NAVES  PROPAGATE  AND  CALCULATE  VERTICAL  TRANSFER  FUNCTION 
FOR  MIC  AND  VERT  GEOPEONE 

SUBROUTINE  TRANSF(B, Ml ,M2 , FR, DPR IKE, D1 ,DMAG, DMMAG, PHAS, PHAS2) 
IMPLICIT  REAL  *8  (A-H,0-Z) 

COMPLEX  *16  J,L1,L2,C1,C2,C3 

COMPLEX  *16  TEMP4 , TPRES , Ml , M2 , TEMF6 

COMPLEX  *16  MOVRM,T1,T2,P1,P2,P3,P4 

COMPLEX  *16  PRES, NV ELI  ,NVEL2,NVEL3  ,NVEL4,NVEL5,NVEL6 

COMPLEX  *16  TNVEL , GOVRM 

COMPLEX  *16  B( 9) 

REAL  *8  M, KFjKBBAR, KR,L0,L3 
COMMON  /BLK7/  H,C,M 

COMMON  /BLK5/  VO, CO, SO f Cl ,L0,L1 ,C2,L2, S3 . L3 ,C3 
COMMON  /BLK3/  RHOF, GB, KBBAR, KR , RHOS, OMEGA, KF 
J  z  ( 0 . DO , 1 . DO ) 

TWOPI  z  4 . DO*DARSIN ( 1 . DO ) 

DPRIME  IS  PROPAGATION  DEPTH 
D1  IS  PROPAGATION  HEIGHT 


CALCULATE  PRESSURE  AT  ABOVE  GROUND  MIC 

PRES  z  -J*KF*L0*( 1 .D0+E( 9) ) 

CALCULATE  NORMAL  COMPONENTS  OF  PARTICLE  VELOCITY 

NVEL1  z  -J*TWOPI*FR*B( 1) *C1*CDEXP(J*L1 *C1*D PRIME) 
NVEL2  z  J*TWOPI*FR*B( 2) *C1*CDEXP(-J*L1 *C1 *D PRIME) 
NVEL3  =  -J*TWOPI*FR*B(3) *C2*CDEXP(J*L2*C2*D PRIME) 
NVEL4  z  J*TWOPI*FR*B( 4) *C2*CDEXP(-J*L2*C2*D PRIME) 
NVEL5  z  J*TWOPI*FR*B(5)*S3*CDEXP(J*L3*C3#DPRIME) 
NVEL6  z  J*TWOPI*FR*B(6)*S3*CDEXP(-J*L3fC3*DPRIME) 
TNVEL  z  NVEL1  +  NVEL2  +  NVEL3  +  NVEL4  +  NVEL5  +  NVEL6 
T1  s  (C  -  M1*M) 

T2  =  (C  -  M2*M) 

PI  z  B(1)*L1*T1*CDEXP(J«L1*C1*DPRIME) 

P2  z  B(2)«L1*T1*CDEXP(-J*L1*C1*DPRIME) 

P3  *  B(3) *L2*T2*CDEXP( J*L2*C2*DPRIME) 

P4  z  B( 4) *L2*T2*CDEXP(-J*L2*C2*DPRIME) 

TPRES  z  PI  +  P2  +  P3  +  P4 


TRA00010 
TR AO 0020 
TRA00030 
TR AO 0040 
TRA00050 
TR AO 0060 
TRA00070 
TRA00080 
TRA00090 
TRA00100 
TRA001 10 
TR  AO  01 20 
TRA00130 
TRA00140 
TR  AO  01 50 
TRA00160 
TRA00170 
TRA001 80 
TR  AO  01 90 
TRA00200 
TRA0021 0 
TR AO 0220 
TRA00230 
TR AO 02 40 
TRA00250 
TR AO 0260 
TR AO 0270 
TRA00280 
TR AO 0290 
TRA00300 
TRA00310 
TRA00320 
TRA00330 
TRA00340 
TRA00350 
TRA00360 
TRA00370 
TRA00380 
TRA00390 
TRA00400 
TRA0041 0 
TRA00420 
TRA00430 
TRA00440 


V’.'V 


W 


J  V  ’J'  •  ■ 


V 


WRITE ( 6,13) PRES , TPRES 
F0RMAT(1I,4(E8.2,2I)) 

MOVRM  =  TPRES/ PRES 
GOVRM  a  TNVEL/PRES 
TEMP4  =  DCONJG(  GOVRM) 

TEMP6  =  DCONJG( MOVRM) 

TEMPS  =  GOVRM«TEMP4 
TEMP7  a  MOVRM«TEMF6 
DMAG  a  DSQRT(TEMPS) 

DMMAG  a  DSQRT(TEMP7) 

TEMPI  a  DIMAG (GOVRM) 

TEMP8  a  DIMAG (MOVRM) 

TEMP2  a  DREAL( GOVRM) 

TEMP9  =  DREAL (MOVRM) 

ARG  a  TEMP1/TEMP2 

ARG2  a  TEMP8/TEMP9 

PH AS  a  DATAN(ARG)«360.D0/TW0PI 

PHAS2  a  DA TAN ( ARG2 ) *360 . DQ/TWOPI 

RETURN 

END 


TRA00450 
TRA00460 
TRA00470 
TRA00480 
TRA00490 
TRA00500 
TRA0051 0 
TRA00520 
TRA00530 
TRA00540 
TRA00550 
TRA00560 
TRA00570 
TRA00580 
TRA00590 
TRA00600 
TRA00610 
TRA00620 
TRA00630 
TRA006 40 


APPENDIX  A. 5 


The  subroutine  LTRANF. FORTRAN  computes  the  vertical  seismic  transfer 
function  in  the  clay  by  taking  the  ratio  of  the  normal  seismic  response  to 
the  acoustic  pressure  at  the  surface.  LTRANF  then  returns  the  magnitude 
(DMAG)  and  phase  (PHAS)  of  the  transfer  function.  Input  arguments  are  all 
defined  in  the  main  routine  AMPL. FORTRAN. 


LET  WAVES  PROPAGATE  AND  CALCULATE  VERTICAL  TRANSFER  FUNCTION 
FOR  MIC  AND  VERT  GEOPHONE  IN  CLAY 

SUB  ROU TINE  LTRAN F( B , LO , RHOF , KF , FR , DPR IME , DMAG , PHAS ) 

IMPLICIT  REAL  *8  (A-H,0-Z) 

COMPLEX  *16  J 

COMPLEX  *16  TEMP4 , NVEL 1 , NVEL2 
COMPLEX  *16  TNVEL , GOVRM, PRES 
COMPLEX  *16  B(9),C4,C5 
REAL  *8  L0,KF,L4,L5 
COMMON  /BLK8/  C4,S4, C5 ,L4,L5 
J  =  (O.DO.I.DO) 

TWOPI  =  4.D0*DARSIN( 1 .DO) 

DPRIME  IS  PROPAGATION  DEPTH 
D1  IS  PROPAGATION  HEIGHT 


CALCULATE  PRESSURE  AT  ABOVE  GROUND  MIC 

PRES  =  -- J*KF*L0  *  ( 1 .  D0+B(  9 ) ) 

CALCULATE  NORMAL  COMPONENTS  OF  PARTICLE  VELOCITY 

NVEL1  =  -J*TW0PI*FR*B(8)*S4*CDEXP(J«L4*C4*DPRIME) 
NVEL2  =  -J*TWOPI*FR*B(7)*C5*CDEXP(J*L5*C5*DPRIME) 
TNVEL  =  NVEL1  +  NVEL  2 
GOVRM  =  TNVEL/ PRES 
TEMP4  =  DCONJG( GOVRM) 

TEMP5  =  GOVRM *TEMP4 
DMAG  =  DSQ  RT(TEMP5) 

TEMPI  =  DIMAG (GOVRM) 

TEMP2  =  DREAL( GOVRM) 

ARG  =  TEMPI /TEMP2 

PHAS  =  DATAN(ARG)»360.D0/TW0PI 

RETURN 

END 


LTR00010 
LTR00020 
LTR00030 
LTR00040 
LTR00050 
LTR00060 
LTR00070 
LTR00080 
LTR00090 
LTROOIOO 
LTR001 10 
LTR00120 
LTR00130 
LTR00140 
LTR00150 
LTR00160 
LTR00170 
LTR00180 
LTR00190 
LTR00200 
LTR00210 
LTR00220 
LTR00230 
LTR00240 
LTR00250 
LTR00260 
LTR00270 
LTR00280 
LTR00290 
LTR00300 
LTR00310 
LTR00320 
LTR00330 
LTR00340 
LTR00350 
LTR00360 
LTR00370 
LTR00380 
LTR00390 
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APPENDIX  A. 6 


RTRANS. FORTRAN  returns  the  manitude  (DMAG)  and  phase  (PHAS)  of  the 
radial  seismic  transfer  function  for  the  porous  sand.  The  routine  is  the 
same  as  TRANSF  except  that  the  radial  component  is  taken  as  opposed  to  the 
normal  component.  The  arguments  in  the  call  statement  are  defined  in  the 
calling  routine,  AMPL. FORTRAN. 


LET  WAVES  PROPAGATE  AND  CALCULATE  RADIAL  TRANSFER  FUNCTION 

SUBROUTINE  RTRANS(B, FR,DPRIME,D1 , DMAG, PHAS) 

IMPLICIT  REAL  *8  (A-H.O-Z) 

COMPLEX  *16  J,L1,L2,C1,C2,C3,S1  ,S2 

COMPLEX  *16  TEMP4 yTPRES,M1 ,M2,TEMP6 

COMPLEX  *16  PRES,jIVEL1  ,NVEL2,NVEL3,NVEL4,NVEL5,NVEL6 

COMPLEX  *16  TNVEL, GOVRM 

COMPLEX  *16  B(9) 

REAL  *8  M,EF,KBBAR,KR,L0,L3 
COMMON  /BLK8/  S1,S2 
COMMON  /BLK7/  H,C,M 

COMMON  /BLK5/  VO.CO.SO ,C1 ,L0 ,L1 ,C2,L2 ,S3 ,L3 ,C3 
COMMON  /BLE3/  RHOF, GB, KBBAR, KR, RHOS, OMEGA, KF 
J  =  ( 0 . DO , 1 . DO ) 

TWOPI  a  4.D0*DARSIN( 1.D0) 

DPRIME  IS  PROPAGATION  DEPTH 
D1  IS  PROPAGATION  HEIGHT 


CALCULATE  PRESSURE  AT  ABOVE  GROUND  MIC 

PRES  a  - J *KF*L0 *  C 1 . D0+B( 9 ) ) 

CALCULATE  RADIAL  COMPONENTS  OF  PARTICLE  VELOCITY 

NVEL1  a  -J*TWOPI*FR*B( 1)*S1  *CDEXP( J*L1 *C1*D PRIME) 
NVEL2  a  J*TWOPI*FR*B( 2)*S1  *CDEXP(-J*L1 *C1*DPRIME) 
NVEL3  a  -J*TWOPI*FR*B( 3) *S2*CDEXP(J*L2*C2*D PRIME) 
NVEL4  a  J*TW0PI*FR*B( 4)*S2*CDEXP(-J*L2*C2*DPRIME) 
NVEL5  =  J*TWOPI*FR*B(5)*C3*CDEXP(J*L3*C3*DPRIME) 
NVEL6  a  J*TWOPI*FR*B(6)*C3*CDEXP(-J*L3*C3*DPRIME) 
TNVEL  =  NVEL1  +  NVEL2  +  NVEL3  +  NVEL4  +  NVEL5  +  NVEL6 

GOVRM  a  TNVEL /PRES 

TEMP4  a  DC0NJG( GOVRM) 

TEMP5  a  G0VRM*TEMP4 

DMAG  a  DSQ  RT(TEMP5 ) 

TEMPI  a  DIMAG  (GOVRM) 

TEMP2  a  DREALC GOVRM) 

ARG  a  TEMPI /TEMP2 

PHAS  a  DATAN(ARG)*360.D0/TW0PI 

RETURN 

END 


RTR0001 0 
RTR00020 
RTR00030 
RTR00040 
RTR00050 
RTR00060 
RTR00070 
RTR00080 
RTR00090 
RTR00100 
RTR001 10 
RTR00120 
RTR00130 
RTR00140 
RTR00150 
RTR00160 
RTR00170 
RTR00180 
RTR00190 
RTR00200 
RTR00210 
RTR00220 
RTR00230 
RTR00240 
RTR00250 
RTR00260 
RTR00270 
RTR00280 
RTR00290 
RTR00300 
RTR00310 
RTR00320 
RTR00330 
RTR00340 
RTR00350 
RTR00360 
RTR00370 
RTR00380 
RTR00390 
RTR00400 
RTR0041 0 
RTR00420 
RTR00430 
RTR00440 
RTR00450 
RTR00460 
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This  routine  was  used  to  compute  the  magnitude  and  phase  of  the 
experimental  transfer  function.  Tranff.Bas  program  first  accepts  the  cali¬ 
brations  (typically  10-15  seconds  depending  on  the  length  of  data  record  to 
be  analyzed)  on  channel  0  (always  the  reference  microphone  above  the 
ground).  The  digitized  record  is  always  collected  in  blocks  of  1024 
points.  Due  to  finite  transfer  rates,  a  few  points  could  be  dropped 
between  successive  1024  point  records.  After  the  calibration  tone  is 
digitized,  each  1024  point  record  is  Fourier  transformed.  The  magnitudes 
of  the  transform  at  each  frequency  (512  of  these)  are  added  together  and 
these  are  added  to  the  magnitude  determined  from  each  successive  1024  point 
data  record.  The  total  sura  is  then  divided  by  the  number  of  data  records  - 
this  is  referred  to  as  the  calibration  level  for  channel  0.  The  tape  is 
then  moved  to  the  calibration  tone  on  the  channel  to  be  analyzed  and  this 
process  is  repeated  to  give  the  calibration  level  for  channel  1. 

Next  the  tape  is  moved  to  the  data  for  the  first  third  octave  band 
(31.5  Hz).  At  this  point,  the  data  from  the  reference  microphone  and 
channel  of  interest  are  simultaneously  digitized.  Actually  channel  0  is 
sampled  first;  channel  one  is  sampled  a  few  hundred  microseconds  later. 
This  small  delay  gives  a  small  phase  difference  which  must  be  corrected  for 
later.  Once  the  data  is  all  digitized  and  stored  in  memory,  each  1024 
point  data  record  from  channel  0  is  Fourier  transformed.  The  simultaneous 
1024  point  record  from  channel  1  is  next  Fourier  transformed  and  the  trans¬ 
fer  function  for  these  1024  points  is  computed  giving  the  transfer  function 
at  512  frequencies.  Only  those  frequencies  within  the  third  octave  band 
being  considered  are  kept.  The  next  pair  of  1024  point  records  is  trans¬ 
formed  and  the  transfer  function  computed  is  added  to  the  first.  This 
process  is  continued  for  all  the  data  records  for  that  third  octave  band. 
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Next,  the  tape  is  advanced  to  the  next  third  octave  band  and  the  process  is 
repeated. 


The  raw  transfer  function  at  each  frequency  is  divided  by  the  calibra-  j 

tion  levels  to  get  absolute  values  for  the  real  and  imaginary  parts  of  the 
transfer  function.  Finally,  the  real  and  imaginary  parts  are  converted  to 
a  magnitude  and  phase  at  each  frequency;  the  results  are  plotted  on  an  X-Y 
plotter  (program  Plotff.Bas,  Appendix  A. 8)  and  stored  on  disks. 

The  plotted  transfer  function  was,  in  many  cases,  quite  noisy  (many 
sharp  spikes).  For  comparison  to  theory,  sometimes  smoothed  this  data.  A  j 

seven  point  averaging  routine  was  used,  twice,  to  smooth  the  data  which  was 
then  replotted.  The  smoothing  routine  (Program  Smooth. Bas,  Appendix  A. 9) 
computed : 

t 

Y(0)  -  {[Y(-3)+Y(3)]+4[Y(-2)+Y(2)]+ll[Y(-l)+Y(l)]+14Y(0)}/42  | 

i 

\ 

\ 

which  replaced  Y(0).  This  procedure  was  executed  twice  to  result  in  graphs 
referred  to  later  as  smoothed  curves. 

As  noted  earlier,  samples  from  the  two  channels  are  actually  taken 
sequentially  which  means  that  there  is  a  small  time  delay  between  samples 
on  the  two  channels.  This  delay  was  determined  experimentally.  An  oscil¬ 
lator  generating  a  sine  wave  at  one  of  the  1/3  octave  band  center 
frequencies  was  coupled  directly  into  channel  0  of  the  computer  and  the 
same  signal  was  put  through  an  RC  filter  (time  constant  .002  sec)  into 
channel  1.  Through  such  a  circuit,  the  amplitude  should  vary  between  0  (at 
0  Hz)  and  1  (well  above  80  Hz,  where  a>  *  1/RC).  The  phase  should  decrease 
from  90°  at  0  Hz  to  0°  well  above  80  Hz.  When  the  amplitude  of  channel  1 
is  0.707  times  the  amplitude  of  channel  0,  the  phase  difference  should  be 
+45°.  Experimentally,  the  phase  was  found  to  be  slightly  larger  than  45° 
requiring  that  we  introduce  a  phase  correction  of  .06  x  frequency.  Note 
that  at  a  frequency  of  6  kHz,  the  phase  corrections  would  be  360°  implying 
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a  channel  Co  channel  delay  of  167  sec  which  is  close  Co  Che  manufacCurer's 


(DEC)  specficaCions.  This  correccion  appears  on  line  5297  of  program 
Tranf f . 


20  OPEN  "VM: FILE1.DAT"  AS  FILE  # 1 
30  OPEN  "VM: FILE2.DAT"  AS  FILE  #2 
40  OPEN  "VM: FILE3.DAT"  AS  FILE  #3 
45  OPEN  "VM: FILE4.DAT"  AS  FILE  #4 
47  OPEN  "LP:"  FOR  OUTPUT  AS  FILE  #7 
50  DIM  #1,A1Z(33768) 

60  DIM  #2,D1(1024) 

70  DIM  #3,F1(1024) 

75  DIM  #4,0(14) 

78  DIM  CZ(1024) 

80  COMMON  AZ(2048) 

85  VZ-0 

100  PRINT  "NUMBER  OF  SECONDS  TO  BE  DIGITIZED?"; 

110  INPUT  S3Z 
120  S4-S3Z-1 

130  PRINT  "LISTEN  FOR  CAL.  TONE  ON  CH."VZ",  PUSH  RETURN"; 

140  INPUT  Z$ 

160  FOR  IZ-0  TO  S4 

170  AIN(,AZ(), 1024, 1/1000, VZ.l) 

180  BL0CK_M0VE( AZ ( 0 ) ,A1Z( IZ*1024 ) , 1024 ) 

190  NEXT  IZ 

210  PRINT  "STOP  TAPE" 

230  S2-0 

240  FOR  IZ-1  TO  512 

250  Dl(IZ)-0 

260  F1(IZ)“0 

270  NEXT  IZ 

280  FOR  IZ=0  TO  S3Z-1 

290  N1“1024*I%  \  N2-N1+1023  \  S9Z»1 

300  GOSUB  7370 

310  FFT( ,1024 ,AZ( ) ,CZ( ) ,SZ) 

320  FOR  JZ-1  TO  512 

330  D1(JZ)-D1(JZ)+SQR((AZ(JZ)*2“SZ)*2-KCZ(JZ)*2*SZ)“2) 

350  NEXT  JZ 
360  NEXT  IZ 
370  FOR  IZ*1  TO  512 
380  S2-S2+D1(IZ) 

390  NEXT  IZ 

400  S2-S2/S3Z 

410  IF  VZ-1  THEN  460 

420  VZ*VZ+ 1 

430  S1-S2 

440  GO  TO  130 

460  JZ-0 

470  PRINT  SI ,S2 

570  0(0)* 20  \  0(l)-25  \  0(2)-31.5  \  0(3)-40  \  0(4)-50  \  0(5)-63  \  0(6)-80 
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670  0(7 >-100  \  0(8)-125  \  0(9)»160  \  0(10)-200  \  0(11)-250  \  0(12)-315 

870  F0-1Q00/1024 

880  FOR  IX- 1  TO  1024 

900  Dl(lX)-0 

910  FI ( IX )-0 

920  NEXT  IX 

930  S4-S3X-1 

970  FOR  LX-1  TO  12 

1070  CLOSE  #1 

1170  OPEN  "VM: FILE 1 .DAT"  AS  FILE  #1 
1270  Q8X-53 

1370  F1-(2~.1667)*0(LX-1) 

1470  F2-(2~.1667)*0(LX) 

1570  I1X-F1/F0 
1670  I2X-F2/F0 

1770  PRINT  "START  TAPE  ";0<LX);M  HZ  AND  PRESS  RETURN."; 

1870  INPUT  D9$ 


SCHMITT( 1,8070) 

IF  Q8XO-999  THEN  2070 

SCHMITT 

FOR  IX-0  TO  S4 

AIN(,AX(), 2048, 1/1000 ,0,2) 

BLOCK_MOVE( AX(0 ) ,  A1X( IX*2048 ) ,2048 ) 

NEXT  IX 

PRINT  "STOP  TAPE  1 !!!!!!!!! !!!  " 

FOR  IX-0  TO  S3X-1 

N1-2*1024*IX  \  N2-N1+2046  \  S9X-2 

GOSUB  7370 

FFT( ,1024 ,AX( ) ,CX( ) ,SX) 

FOR  MX-I1X  TO  I2X 
JX-MX+5I2 

D1(MX)-(AX(MX)*2~SZ)/S1 
D1(JX)-(CX(MZ)*2~SX)/S1 
NEXT  MX 

N1-2*1024*IX+1  \  N2-N1+2046  \  S9X-2 
GOSUB  7370 

FFT( ,1024 ,AX( ) ,CX( ) ,SX) 

FOR  MX-I1X  TO  12 X 
JX-MX+512 

M-Dl  (MX)',2+D1  C  JX)',2 

R-(AX(MX)*D1(MX)+CX(MX)*D1(JX))*2“SX/S2 
I2-(D1(MX)*CX(MX)-A%(M%)*D1(JX))*2~SX/S2 
FI  (MX)-Fl  (MX  )+SQR<  (R~2+I2*2  )/M“2  ) 

FI ( JX )-Fl ( JX )+ATN( 12 /R) 

NEXT  MX 
NEXT  IX 

FOR  IZ-I1X  TO  12 X 
MX-IX+512 
D1(IX)-IX*F0 
F1(IX)-F1(IX)/S3X 
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5297  F1(M2)-180*((F1(M2)/S32)/3.14159)-D1(I2)*60/1000 
5310  NEXT  12 
5370  NEXT  L2 

5870  II 2*  (2'“  .1667 )*0(0)/F0 
5970  I22*(2~.1667 )*0(13)/F0 
6070  FOR  12-112  TO  122 
6170  M2* 12+5 12 

6780  PRINT  #7  ,"FREQ*"D1(I2)  ,"MAG-"F1 ( 12 )  , "PHASE* "FI (M2 ) 
6970  NEXT  12 
7229  CLOSE 

7250  CHAIN  "PLOTFF.BAS" 

7270  STOP 
7370  J2-0 

7380  FOR  I32*N1  TO  N2  STEP  S92 
7470  J2-J2+1 
7570  a(J2)*0 
7670  A2(J2)«A12(I32) 

7770  NEXT  132 
7870  RETURN 
7970  END 
8070  Q82-- 999 
8170  RETURN 
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PLOTFF. BAS  was  written  to  plot  the  output  of  TRANFF. BAS  on  an  HP  7470A 
digital  plotter. 


60  OPEN  "VM: FILE2.DAT"  AS  FILE  #2 
70  OPEN  "VM: FILE3.DAT"  AS  FILE  # 3 
80  DIM  #2 ,DL ( 1024 ) 

90  DIM  #3, FI (1024) 

100  L9-2  \  0-1  \  01-1  \  Q-l 
200  DISPLAYjCLEAR 

1400  PRINT  "  input  the  size  of  plot  (1)  full  page,  and  (2)  quarter  page  ? 
1500  INPUT  AZ 
1600  X9-0  \  X8-0 
3400  FOR  1-1  TO  2 

3450  PRINT  "PLEASE  PUT  PAPER  IN  H.P.  PLOTTER  AND  PRESS  RETURN";  \  INPUT  E7$ 
3475  Y9-0  \  Y8-0  \  L9-2  \  0-1  \  01-1  \  Q-l 
3500  FORJ-l  TO  362 

3600  IF  1-1  THEN  IF  Y9<F1(J)  THEN  Y9-F1(J) 

3700  IF  1-1  THEN  IF  Y8>F1(J)  THEN  Y8-F1(J) 

3800  IF  1-2  THEN  IF  Y9<F1(J+512)  THEN  Y9-F1(J+512) 

3900  IF  1-2  THEN  IF  Y8>F1(J+512)  THEN  Y8-F1(J+512) 

3950  IF  1-2  THEN  4500 

4000  IF  X9<D1( J)  THEN  X9-D1(J) 

4100  IF  X8>D1(J)  THEN  X8-D1(J) 

4500  NEXT  J 

4600  IF  Y9<3  THEN  01-10000 
4700  IF  X9<3  THEN  0-10000 
5100  PAUSE( 3 ) 

5200  IF  Y9*01<0  THEN  PI— 1 
5300  IF  Y9*01>-0  THEN  Pl-1 
5400  Y2»INT( L0G10(ABS( Y9*01 ) ) ) 

5500  IF  Y2<0  THEN  Y2-Y2-1 
5600  Y3-Y9*01/10“Y2 
5700  IF  ABS(Y3)-1  THEN  Sl-1 
5800  IF  ABS(Y3 )> 1  THEN  Sl-2 
5900  IF  ABS(Y3)>2  THEN  Sl-2. 5 
6000  IF  ABS(Y3 )> 2.5  THEN  Sl-5 
6100  IF  ABS(Y3 )>5  THEN  Sl-10 
6200  U1-(S1*10~Y2)*P1 

6300  IF  ABS(Y8*01)<(1.00000E-03*(Y9*01))  THEN  7600 
6400  IF  Y8*01<0  THEN  P2— 1 
6500  IF  Y8*01>-0  THEN  P2-1 
6600  Y4-INT(L0G10(ABS( Y8*01 ))) 

6700  IF  Y4<0  THEN  Y4-Y4-1 
6800  Y5-Y8*01/10~Y4 
6900  IF  ABS(Y5 )-l  THEN  S2-1 
7000  IF  ABS( Y5 )> 1  THEN  S2-2 
7100  IF  ABS(Y5)>2  THEN  S2-2.5 
7200  IF  ABS(Y5)>2.5  THEN  S2-5 
7300  IF  ABS(Y5)>5  THEN  S2-10 
7400  Ll«( S2*10“Y4 )*P2 
7500  GO  TO  8300 
7600  Ll-0 
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8200  PAUSE(3) 

8300  IF  AX-1  THEN  X6-9200  \  Y6-7560 
8400  IF  AX-2  THEN  X6-5150  \  Y6-4180 
8410  IF  ABS(U1)> 30000  THEN  Q-100 
8420  IF  ABS(L1  )>30000  THEN  Q-100 
8500  C0UT( ,"IN;SP1 ,1 ) 

8600  M$— "IP  1100,1 600, "+STR$(X6)+","+STR$(Y6)+";" 

8700  C0UT( ,M$ , ,1 ) 

8800  M$-"SC  "+STR$ ( X8*0 )+" ,"+STR$ ( X9*0 )+" , "+STR$ ( LI /Q)+" , "+STR$ ( U1 /Q )+" ; " 
8900  C0UT( ,M$ , ,1 ) 

9000  M$-"PA  "+STR$ ( X8*0 )+" , "+STR$ ( LI /Q )+" , PD , "+STR$ ( X9*0 )+" , "+STR$ ( LI /Q )+" 
9 100  M$-M$+STR$  (  X9*0  )+" ,  "+STR$  (  U1  /Q )+" ,  "+STR$  (  X8*0  )+" ,  "+STR$  (  U1  /Q )+" , 
"+STR$ (X8*0)+" , "+STR$ (Ll/Q)+" , PU ; " 

9200  C0UT(,M$,,1) 

9300  PAUSE( 2 ) 


9900  C0UT( , "SR  1,2;",,1) 

10000  Tl-(ABS(X9*0)+ABS(X8*0))/5 
10100  T2-(ABS(Ll/Q)+ABS(Ul/Q))/5 
10200  C-X8*0 
10300  FOR  J-l  TO  6 

10400  M$-"PA  "+STR$ ( C )+" , "+STR$ ( LI /Q )+" ;XT ; " 

10500  C0UT( ,M$ , ,1 ) 

10600  M$-"PR  "+STR$ (-(Tl*. 15 ))+","+STR$(-(T2*.2 ))+";” 

10700  C0UT( ,M$ , ,1 ) 

10800  M$-"LB"+STR$(C/0)+CHR$(3) 

10900  C0UT( ,M$ , ,1 ) 

10950  PAUSE( 2 ) 

11000  C-C+Tl 
11100  NEXT  J 
11200  PAUSE(4 ) 

11300  C-Ll/Q 
11400  FOR  J-l  TO  6 

11500  M$-"PA  "+STR$ ( X8*0 )+" , "+  STR$ ( C )+" ; YT ; " 

11600  C0UT( , M$ , , 1 ) 

11700  M$-"PR  "+STR$(-(T1*.4 ) )+" ,0 ;LB"+STR$ ( (C/01 )*Q)+CHR$(3) 

11800  C0UT( , M$ , , 1 ) 

11900  C-C+T2 
12000  NEXT  J 
12050  PAUSE( 5 ) 

12600  M$-"PA  "+STR$(X8*0)+",0 ;PR  "+STR$(-(T1*.4 ) )+" ,0 ;LB0  "+CHR$(3) 
12700  C0UT( ,MS , ,1 ) 

12800  M$-"PA  "+STR$(X8*0)+",0,PD,"+STR$(X9*0)+",0,PU;" 

12900  C0UT( ,M$ , ,1 ) 

13000  PAUSE(5) 

13100  IF  1-1  THEN  A2$ -"MAGNITUDE" 

13200  IF  1-2  THEN  A2$-"PHASE" 

13300  A1$-"FR£QUENCY" 

14000  Z1-LEN(A1$) 

14100  Z2-LEN( A2$ ) 

14200  M$-"PA  "+STR$ ( ( ( X8*0 )+( X9*0 ) ) /2 )+" , "+STR$ ( LI /Q)+" ; " 


14300  C0UT( ,M$ , ,1 ) 

14400  M$-"PR  0,"+STR$(-(T2*. 4 ))+";" 

14500  C0UT( ,M$ , ,1 ) 

14600  PAUSE(l) 

14700  M$-"PR  "+STR$(-(Tl*.05)*(Zl/2))+",0;LB"+Al$+CHR$(3) 

14800  C0UT( ,M$ , ,1 ) 

14900  PAUSE( 2 ) 

15000  M$-"PA  "+STR$(X8*0)+","+STR$ ( ( (LI /Q)+(U1  /Q) )/2 )+"; " 

15100  C0UX( ,M$ , ,1 ) 

15200  M$-"PR  "+STR$ (  -(  T1  * .  6 ) )+"  ,0 ; " 

15300  C0UT( ,M$ , ,1 ) 

15400  M$«"PR  0 ,"+STR$(-(T2*.06 )*(Z2/2 )  )+" ;DR  0 ,1 ;LB"+A2$+CHR$(3) 

15500  C0UT(,M$,,1) 

15550  C0UT(,"SP2;",,1) 

16200  PAUSE(20 ) 

17200  M$»"PA  "+STR$(D1(1 )*0)+","+STR$((Fl ((1-1 )*512+1 )*01 )/Q)+";PD;" 

17300  C0UT(,M$,,1) 

17400  FOR  >2  TO  362 

17450  IF  INT(J/100)-J/100  THEN  PAUSE(l) 

17600  M$-"PA  "+STR$ ( D1 ( J )*0 )+" , "+STR$ ( ( FI ( ( ( I— 1 )*5 12 )+J )*01 ) /Q )+" ; " 

17700  C0UT(,M$,,1) 

17800  NEXT  J 

18000  C0UT(,"PU;",,1) 

18100  PAUSE(3) 

18200  NEXT  I 

18210  PRINT  "WOULD  YOU  LIKE  TO  SAVE  THESE  NUMBERS  ON  DISK?  IF  SO  TYPE  A  FILE" 
18220  PRINT  "NAME  IN.  IF  YOU  DON’T  WANT  TO  SAVE  THESE  NUMBER  JUST  PRESS  RETURN 
18230  PRINT  "INPUT  FILE  NAME";  \  INPUT  F4$ 


18240  IF  F4$-""  THEN  18300 

18250  OPEN  "DY1 : "+F4$+".XYZ"  AS  FILE  #1 

18260  FOR  LZ-1%  TO  362% 

18270  PRINT  //I  ,D1  (L%)","F1  (L%)","F1  (L%+5 12%  ) 
18280  NEXT  L% 

18290  CLOSE 
18300  END 
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SMOOTH. BAS  is  a  seven  point  smoothing  routine  used  to  smooth  the 


output  of  TRANFF.BAS  and  is  discussed  in  Chapter  4  of  this  report. 


10  OPEN  "VM: FILE2.DAT"  AS  FILE  #2 
20  OPEN  "VM: FILE3.DAT"  AS  FILE  #3 
30  DIM  #2,D1(1024) 

40  DIM  #3,F1(1024) 

50  DIM  X(512),Y(512),Z(512) 

60  PRINT  "INPUT  FILE  NAME";  \  INPUT  Al$ 

70  OPEN  "DY1 :"+Al$+".XYZ"  AS  FILE  #1 

80  FOR  IZ-1  TO  353 

90  INPUT  #1,X(IZ),Y(IZ),Z(IZ) 

100  NEXT  I Z 
110  NZ-0 

120  FOR  IZ-4  TO  350 

130  Y( IZ)-Y( IZ)*14+1 1*(Y( IZ-1 )+Y( IZ+1 ) )+4*(  Y( IZ-2 )+Y( IZ+2) ) 

140  Y( IZ )«< Y( IZ)-( Y( IZ-3 )+Y( IZ+3 ) ) )/42 

141  IF  Z(IZ)<90  GO  TO  143 

142  Z(  IZ)»90 

143  IF  Z(IZ)>-90  GO  TO  150 

144  Z(IZ)«- 90 

150  Z( IZ )«Z< IZ)*14+11*(Z( IZ-1 )+Z< IZ+1 ) >+4*(Z( IZ-2 )+Z( IZ+2 ) ) 
160  Z( IZ )-(Z( IZ)-(Z( IZ-3 )+Z< IZ+3 ) ) )/42 

170  NEXT  IZ 

171  Y(350)-0 

172  Z(350)-0 
180  NZ-NZ+1 

190  IF  NZ-1  GO  TO  120 

200  BLOCK_MOVE(X(0),D1(0),512) 

210  BL0CK_M0VE( Y(0 ) ,F1 (0 ) ,512 ) 

220  BLOCK_M0VE(Z(0 ) , FI (513), 5 12) 

230  CLOSE 

240  CHAIN  "PLOTFF . BAS" 

250  END 
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Included  here  is  a  larger  and  more  complete  sample  of  the  experimental 
data  collected.  The  purpose  is  twofold;  one,  to  present  a  sample  of  the 
raw  data  before  smoothing  calculations  were  done  and  two,  to  make  available 
samples  of  the  radial  and  transverse  seismic  transfer  function  which  were 
not  included  in  the  body  of  the  report.  The  data  chosen  for  presentation 
were  collected  at  the  UM  test  field  and  have  not  been  scaled.  The  units 
have  been  left  off  the  axes,  however,  they  are  the  same  as  on  earlier 
graphs.  On  each  figure  is  a  legend  which  includes  the  date,  angle  of 
incidence,  depth  of  the  sensor  and  the  axis  of  the  geophone  or  probe. 

Figures  1-4  indicate  the  poor  signal  to  noise  ratio  obtained  as  the 
angle  of  incidence  was  changed  from  20°  to  5°.  The  remaining  eight  figures 
are  the  magnitudes  and  phases  of  the  vertical,  radial  and  transverse  seis¬ 
mic  transfer  function  and  the  probe  transfer  function. 
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APPENDIX  C.2 


Sardis  beach  is  located  in  the  John  W.  Kyle  State  Park,  which  is 
approximately  10  miles  west  of  Oxford,  Mississippi.  The  experiment  was 
carried  out  on  the  lower  lake  beach  immediately  below  the  eastern  most  end 
of  the  dam.  The  experimental  site  is  marked  with  an  X. 
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